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1  Introduction  and  Objectives  of  Phase  I  Research 

The  problem  of  detection,  classification  and  localization  of  multiple  ground  targets,  e.g.  trucks  and  tanks, 
using  arrays  of  unattended  passive  acoustic  sensors  is  complicated  due  to  various  factors.  These  include: 
variability  and  nonstationarity  of  target  acoustic  signatures,  signal  attenuation,  loss  of  coherence  and  fading 
effects  as  a  function  of  range  and  doppler,  effects  of  wind  noise  and  weather/environmental  conditions  on 
the  sensory  data,  competing  clutter  with  similar  spectral  characteristics,  and  lack  of  a  priori  knowledge 
about  the  number  of  the  targets  in  the  battlefield.  In  addition,  the  presence  of  multiple  closely  spaced 
targets  can  significantly  complicate  the  detection/classification,  data  association  and  localization  processes, 
especially  when  the  tracks  of  the  targets  are  allowed  to  cross,  split,  or  remain  close  together.  As  a  result, 
new  schemes  are  needed  to  provide  fast  and  extremely  accurate  bearing  angle  estimation  and  subsequent 
classification  and  localization  of  different  types  of  narrowband  and  wideband  sources  from  their  acoustic 
signatures.  However,  in  developing  bearing  angle  estimation  and  classification  schemes  for  this  particular 
problem  one  has  to  consider  power  consumption  and  size  constraints,  which  in  turn  place  severe  constraints 
on  the  complexity  of  the  DSP  and  communication  systems  for  these  arrays. 

Unattended  passive  acoustic  sensors  [l]-[3]  are  among  the  widely  used  sensors  for  remote  battlefield 
surveillance,  situation  awareness  and  monitoring  applications.  Together  with  dedicated  digital  signal  pro¬ 
cessing  (DSP),  these  small  and  cost  effective  sensors  can  provide  real-time  information  about  different 
types  of  ground  and  airborne  targets.  They  are  rugged  and  reliable  and  can  be  left  in  the  field  for  a  long 
period  of  time  after  deployment.  The  DSP  associated  with  these  sensors  must  be  capable  of  performing 
various  tasks  including  target  detection  and  direction  of  arrival  (DOA)  estimation,  tracking,  classification 
and/or  identification.  Generally,  there  can  be  a  wide  variety  of  target  types  in  battlefield  depending  on 
the  specific  mission,  e.g.  troops,  ground  targets  (trucks,  tanks,  etc),  airborne  targets  (helicopters,  missiles, 
airplanes)  with  signatures  that  overlap  both  spectrally  and  temporally.  In  addition,  optimum  performance 
is  highly  dependent  on  terrain,  weather,  and  background  noise. 

The  work  at  the  Army  Research  Laboratory  (ARL)  [2]- [3]  involved  development  of  array  signal  process¬ 
ing  algorithms  using  small  baseline  acoustic  arrays  to  estimate  DOA’s,  track,  and  classify  moving  targets 
from  their  acoustic  signatures.  Since  conventional  beamforming  schemes  generally  do  not  provide  good 
results  in  these  situations,  owing  to  the  structural  limitations  of  the  array  and  lack  of  spatial  coherence, 
high-resolution  DOA  algorithms  must  be  developed.  Current  efforts  on  incoherent  and  coherent  wideband 
MUSIC  (Multiple  Signal  Classification)  show  some  promise  for  detecting  and  tracking  multiple  ground 
vehicles.  In  [2]  by  exploiting  the  multi-spectral  content  of  the  targets,  better  results  were  reported  using 
wideband  signal  subspace  DOA  algorithms.  Experimental  results  using  a  circular  array  of  six  sensors, 
with  a  diameter  of  8  ft  indicated  the  advantages  of  the  wideband  MUSIC  over  the  delay-sum  algorithm. 
In  [3]  focused  wideband  adaptive  array  processing  algorithms  are  employed  for  the  high-resolution  DOA 
estimation  of  ground  vehicles.  Several  methods  including  steered  covariance  matrix  (STCM)  and  spatial 
smoothing  or  interpolation  using  the  array  manifold  interpolation  are  considered.  In  addition,  this  study 
provides  experimental  analysis  of  incoherent  and  coherent  wideband  MUSIC  algorithms  and  narrowband 
MUSIC  for  a  circular  array  of  acoustic  sensors  in  terms  of  accuracy  and  stability  of  DOA  estimates.  It 
was  shown  that  for  certain  signal-to-noise  ratio  (SNR)  conditions  incoherent  wideband  methods  yield  more 
accurate  DOA  estimates  than  the  coherent  methods  for  highly  peaked  spectra  while  for  sources  with  flat 
spectra  the  coherent  wideband  methods  generate  more  accurate  DOA  estimates.  The  study  are  indicated 
that  the  coherent  wideband  methods  are  much  more  statistically  stable  than  the  incoherent  ones.  The 
price  paid  for  this  improved  stability  is  high  computational  cost. 

As  far  as  we  know,  none  of  the  existing  algorithms  demonstrated  proven  capability  to  detect,  classify 
and  track  multiple  closely  spaced  sources  in  tight  formations  especially  in  presence  of  clutter  and  wind 
noise.  This  is  critical  in  realistic  battlefield  scenarios  where  an  overwhelmingly  large  number  of  targets  may 
be  encountered  under  changing  and  difficult  environmental/terrain  conditions.  In  addition,  most  of  the 
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existing  methods  apply  to  arrays  with  linear  placement  of  sensors  and  uniform  spatial  sampling.  The  goal  of 
this  Phase  I  research  effort  was  to  study  and  develop  new  DOA  estimation  and  target  classification  methods 
that  can  be  used  to  isolate  multiple  sources  and  classify  them  based  upon  the  spectral/temporal  features  in 
their  acoustic  signatures  as  well  as  their  dynamical  behavior.  We  have  carefully  studied  the  geometry  of  the 
arrays  and  the  properties  of  data  sets  and  developed  two  classes  of  algorithms  for  localizing  and  resolving 
multiple  closely  spaced  targets.  The  first  class  of  approaches  is  based  upon  high-resolution  subspace-based 
DOA  estimation  methods  that  can  be  used  for  localizing  and  tracking  the  targets.  The  second  class, 
however,  is  based  upon  source  separation  methods  using  independent  component  analysis  (ICA)  and  its 
variants  for  determining  the  number  of  potential  targets.  The  primary  criterion  was  to  develop  algorithms 
that  can  potentially  be  implemented  in  battlefield  operations.  In  addition,  the  methods  should  not  rely 
on  prior  knowledge  of  the  number  of  targets,  target’s  dynamical  information  and  initial  conditions,  and 
background  interference  and  clutter.  Of  particular  importance  was  to  investigate  how  the  detection  and 
extraction  of  useful  target  information  can  be  done  in  spatially  aliased  data  due  to  sparse  nature  of  the 
arrays  employed  for  this  problem.  Testing  of  the  algorithms  have  been  done  on  the  real  data  acquired  from 
US  Army  TACOM-ARDEC  in  Picatinny  Arsenal,  NJ  as  well  as  on  synthesized  data  containing  multiple 
target  scenarios  for  a  circular  array  with  a  wagon- wheel  pattern  geometry  [2],  The  effectiveness  of  the 
different  schemes  were  demonstrated  and  benchmarked  on  these  multiple  target  scenarios. 

The  organization  of  this  final  report  is  as  follows.  Section  2  provides  a  review  of  all  the  existing 
methods  for  multiple  wideband  source  separation  and  DOA  estimation  problems  using  passive  acoustic 
arrays.  The  barrier  issues  in  developing  efficient  algorithms  for  detection,  classification  and  tracking  of 
multiple  closely  spaced  targets  in  realistic  battlefield  conditions  are  also  introduced.  Section  3  introduces 
two  different  methods  for  generating  realistic  synthesized  multiple  target  cases  with  various  tight  formations 
e.g.  staggered,  abreast  and  single-file.  Section  4  investigates  and  develops  various  efficient  wideband  DOA 
estimation  methods  that  consider  the  geometry  of  the  arrays  for  high-resolution  and  robust  target  bearing 
angle  estimation  and  separation.  In  particular,  we  investigated  how  to  resolve,  with  minimum  ambiguity 
and  high  confidence,  the  DOA’s  associated  with  multiple  closely  spaced  targets  in  spatially  aliased  data 
due  to  sparse  nature  of  the  arrays  employed  for  this  problem.  The  effectiveness  of  the  algorithms  are 
demonstrated  on  both  real  and  synthesized  data.  The  algorithms  are  then  benchmarked  in  terms  of  their 
accuracy,  resolution  and  robustness  in  estimating  DOA’s  and  separating  multiple  closely  spaced  targets.  In 
Section  5  we  present  some  preliminary  results  of  a  dual-size  array  manifold  that  consists  of  a  rectangular 
array  grid  to  be  spaced  much  greater  than  a  half-wavelength,  but  at  each  grid  has  a  circular  array  geometry. 
Simulation  results  on  a  narrowband  case  are  also  presented.  Section  6  studies  and  develops  various  source 
separation  algorithms  that  exploit  ICA  methods  using  higher  order  statistical  properties  of  the  sources. 
Different  ICA  algorithms  are  extended  and  applied  to  this  problem.  Results  on  several  target  scenarios  are 
also  presented.  The  performance  of  all  the  methods  considered  in  this  Phase  I  research  in  terms  of  accuracy, 
robustness  and  computational  efficiency  are  compared  and  the  promising  methods  for  future  investigation 
have  been  identified.  Section  7  proposes  a  new  algorithm  for  combining  and  associating  the  DOA  estimates 
from  multiple  nodes  in  order  to  localize  and  identify  all  the  sources  in  the  formation.  Finally,  Section  8 
gives  the  concluding  remarks  and  important  topics  for  Phase  II  research. 

2  Identification  of  Barrier  Issues  &  Review  of  Existing  Methods 

In  this  section  the  main  barrier  issues  in  the  development  of  robust,  high-resolution  and  real-time  algo¬ 
rithms  for  tracking  and  classifying  multiple  targets  in  realistic  environment  operating  conditions  using 
unattended  passive  acoustic  arrays  are  briefly  presented.  Below  we  provide  a  review  of  some  of  the  existing 
methods  and  present  the  most  critical  issues,  based  on  the  state-of-art  knowledge  in  the  literature,  that  are 
considered  worthwhile  for  further  investigation.  It  is  our  opinion  that  the  issues  indicated  below  present 
some  of  the  most  challenging  directions  to  be  researched  in  the  next  stages  of  this  project  and  in  the  future. 
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2.1  Narrowband  Direction-of- Arrival  (DOA)  Estimation 

The  problem  of  DOA  estimation  from  a  sensor  array  has  received  much  attention  in  recent  years.  This 
problem  has  been  extensively  studied,  especially  for  narrowband  sources,  a  large  body  of  literature  exists. 
The  most  basic  configuration  for  a  uniform  linear  array  (ULA)  consists  of  M  identical  and  omnidirectional 
sensors  that  receive  ( N  <  M)  narrow-band  signals  from  the  sources  with  unknown  DOA’s  #i,  •  •  •  ,  9pj.  If 
the  signals  are  transformed  to  baseband,  the  received  signal  vector  for  the  nth  source  can  be  written  as 

y(t)  =  a(9n)s(t)  +  v(t),  t  =  l,---,L 

where  a(9n )  is  the  steering  vector  given  by 

a(9n)  =  [1  e~iu)cT2n,  •  •  •  ,  e~^™n}T 

where  Tkn  is  the  time  delay  associated  with  the  kth  sensor  given  by  Tkn  =  with  c  being  the  velocity 

of  sound  propagation  in  air,  9n  is  the  DOA  of  the  nth  source,  and  dk  is  the  distance  between  the  kth 
sensor  and  the  1st  sensor  (reference).  Note  that  it  is  assumed  that  dx  =  0.  Furthermore,  for  ULA  we  have 
dk  =  (k  —  1  )d.  For  N  sources,  the  above  equation  becomes 

Y(t)  =  A(6)S(t)+r(t) 

where  Y(t),  S(t)  and  X (£)  are  the  vectors  of  received  signals,  source  signals  and  additive  noise,  respectively 
and  A(6)  is  the  composite  steering  matrix  of  size  M  x  N  given  by 

A(6)  =  [a(0i),---  ,a(0jv)] 

where  9  =  {6X,---  ,  0jv]T  is  the  DOA  vector.  If  the  additive  noise  is  assumed  to  be  zero-mean  spatially  and 
temporally  white  Gaussian  process  with  the  unknown  diagonal  correlation  matrix  Q  =  d,iag{o\,  ■  •  •  ,  cr2M}, 
the  correlation  matrix  of  the  received  signal  vector,  Y (£) ,  can  be  written  as 

Ry  =ARsAh  +  Q 

where  Rs  =  E{S(t)SH(t)}  is  the  source  correlation  matrix  and  H  denotes  the  Hermitian  transpose. 

The  correlation  matrix  Ry  can  be  computed  from  the  baseband  data  using  Ry  —  ^  J2t=i  ^ 
where  L  is  the  number  of  snapshots  in  one  measurement.  Assuming  that  the  sources  are  uncorrelated 
(incoherent)  and  the  Rank(ARsAH )  =  N,  then  the  eigenvalues  of  matrix  Ry  can  be  arranged  in  de¬ 
creasing  order  i.e.  Ai  >  A2  >  •••A n  »  Av+i  >  •••  >  A m-  The  orthonormal  eigenvectors  associated 
with  {Ai,  A2,  •  •  •  XN},  and  {A^+i,  A^+2,  •  •  •  AM}  form  the  signal  and  noise  subspaces.  If  we  denote  these 
subspaces  by  Gs  =  {ex,  e2,  •  •  •  eN},  and  Gn  =  {eN+i,eN+2,  •  •  •  eM},  respectively,  where  e,  is  the  eigenvector 
associated  with  eigenvalue  A i,  then  we  can  easily  show  that  AHGn  =  0  i.e.  the  columns  of  Gn  belong  to 
the  null  space  of  AH .  This  leads  to  the  important  result  that  the  DOA’s,  0n’s  are  the  only  solutions  of  the 
equation 

aH(9n)GnG”a(9n)  =  0 

Using  the  MUSIC  algorithm,  DOA’s  can  be  determined  as  the  locations  of  the  N  highest  peaks  of  the 
function,  for  6n  e  [— 7r,7r],  of  the  function 

P{9n)  =  atf(0n)GnG£a(0n) 

To  avoid  ambiguity  in  estimating  the  DOA’s  in  MUSIC,  the  Nyquist  criterion  is  imposed,  which  gives 
\ucTmax\  <  7T.  Clearly,  this  limits  the  use  of  MUSIC  for  wideband  sources  using  a  sparse  array  as  spatial 
aliasing  will  be  inevitable  in  such  cases. 
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Given  that  the  elements  of  a  sensor  array  are  commonly  organized  in  a  very  regular  and  structured 
geometry  such  as  linear,  circular,  and  rectangular  shapes,  many  researchers  have  exploited  these  geometrical 
structures  to  develop  computationally  efficient  algorithms  for  the  DOA  estimation.  For  instance,  Root- 
MUSIC,  MODE,  and  IQML  [4]  are  among  some  of  the  well-known  techniques  that  take  advantage  of  a 
uniform  structure  in  the  array  field.  ESPRIT  is  an  alternative  algorithm  that  requires  an  array  consisting  of 
two  identical,  translated  sub-arrays  to  determine  the  DOA  estimates.  Of  particular  interest  to  this  project 
is  the  so-called  MI-ESPRIT  [5]  algorithm  that  attempts  to  generalize  the  standard  ESPRIT  algorithm 
to  arrays  that  are  composed  of  multiple  identical  sub-arrays  arranged  by  translations  of  sub-arrays  with 
arbitrary  but  known  spacings.  More  details  are  provided  in  the  subsequent  subsection. 

2.2  Wideband  DOA  Estimation 

For  wideband  sources,  however,  the  literature  is  somewhat  limited.  Some  of  the  earlier  work  in  this 
domain  are  due  to  Wax  et  al.  [6],  Wang  and  Kaveh  [7],  Buckley  and  Griffiths  [8],  and  Su  and  Morf  [9], 
Both  incoherent  aggregation  technique  suggested  in  [6]  and  coherent  signal  subspace  technique  proposed 
in  (Wang  &  Kaveh;  [7])  represent  in  some  way  transition  of  narrowband  processing  techniques  to  the 
wideband  cases.  Furthermore,  Bienvenu  [10]  and  Buckley  and  Griffiths  [8]  have  shown  how  to  divide  the 
spatio-temporal  observation  space  into  signal  and  noise  subspaces  and  have  suggested  search  functions  to 
estimate  the  DOA’s.  Grenier  [11]  has  generalized  this  notion  to  a  so-called  filter-spatial  domain.  Su  and 
Morf  [9]  have  suggested  a  modal  decomposition  of  the  signal  subspace  for  estimating  the  DOA’s. 

The  problem  of  wideband  DOA  estimation  based  on  maximum  likelihood  principle  is  also  investigated 
in  [12].  The  authors  in  [12]  have  proposed  several  search  functions  to  obtain  the  ML  estimate  of  DOA 
for  wideband  sources.  Due  to  its  extensive  and  time-consuming  nature  this  approach  is  computationally 
expensive  and  impractical  for  real-time  applications.  The  approach  that  is  proposed  by  Bresler  and  Ma- 
covski  [13]  for  the  narrowband  case,  and  a  uniform  linear  array  (ULA),  formulates  the  problem  as  an 
iterative  method  known  as  the  “iterative  quadratic  maximum  likelihood”  (IQML)  for  obtaining  the  ML 
solution.  The  IQML  method  involves  the  creation  of  a  Toeplitz  matrix  that  spans  the  subspace  orthogonal 
to  the  signal  subspace  and  can  be  parameterized  by  the  coefficients  of  a  prediction  polynomial,  whose 
roots  yield  the  DOA  estimates.  In  this  case,  the  use  of  a  parameterized  model  is  equivalent  to  the  prob¬ 
lem  of  estimating  the  coefficients  of  this  model  from  which  DOA’s  can  easily  be  calculated.  In  [14]  the 
authors  have  developed,  based  on  the  ideas  available  for  the  parameter  estimation  of  multiple  narrowband 
sources  in  WGN,  pseudo-maximum  likelihood  parameter  estimators  for  the  wideband  sources  by  adapting 
the  harmonic  model  of  Kay  et  al.  [15].  They  have  shown  that  the  concept  of  an  ARMA  model  for  the 
observed  data  vector  used  in  the  narrowband  case  [13]  can  be  generalized  to  model  the  spatio-temporal 
samples.  As  in  the  narrowband  case,  it  is  shown  that  the  coefficients  of  the  resulting  ARMA  model  can 
be  used  to  represent  the  null  subspace  of  the  wideband  array  steering  matrix.  This  has  resulted  in  the 
formulation  of  an  IQML-like  procedure  for  the  parameter  estimation  of  wideband  sources  that  is  similar  to 
that  for  the  narrowband  case  as  suggested  in  [13],  The  proposed  approach  is  applicable  to  both  coherent 
and  incoherent  sources.  The  performance  of  the  proposed  method  was  also  studied  on  simulated  cases. 

2.3  DOA  Estimation  for  Non-uniform  Arrays 

The  problem  of  DOA  estimation  for  non-uniform  linear  array  (NLA)  has  also  been  investigated  in  [16-22]. 
Nrvn -uniformly  spaced  (or  sparse)  arrays  are  still  the  subject  of  active  investigation  for  both  theoretical 
and  practical  reasons.  Theoretically,  it  has  been  known  for  many  years  that  non-uniform  spacing  may 
lead  to  significantly  improved  DOA  estimation  performance  [23],  Research  on  the  design  of  “optimal” 
non-uniform  geometry  is  still  in  progress.  In  addition  to  DOA  estimation  accuracy  issues,  the  problems 
with  identifiability  and  ambiguity  are  also  under  active  investigation  [24] .  A  crucial  requirement  for  DOA 
estimation  by  subspace-based  methods  is  the  assumption  of  linear  independence  among  the  steering  vectors. 
In  [25]  Schmidt  has  classified  manifold  ambiguities  according  to  their  “rank”  based  on  the  algebraic  linear 
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dependence  of  these  vectors.  Consequently,  most  DOA  estimation  techniques  based  on  subspace  approach 
have  avoided  such  manifold  ambiguity  scenarios.  For  example,  the  study  of  the  unique  localization  of 
multiple  sources  by  passive  sensor  arrays  (identifiability)  conducted  by  Wax  and  Ziskind  [26]  has  been 
based  on  the  assumption  that  “any  subset  of  distinct  steering  vectors  from  the  array  manifold  is  linearly 
independent.”  In  this  paper,  identifiability  conditions  have  been  defined  in  terms  of  the  number  of  sensors 
and  the  rank  of  the  inter-source  correlation  matrix,  under  this  basic  assumption.  As  a  consequence  of 
Schmidt’s  definition,  the  ambiguity  problem  has  usually  been  associated  with  manifold  ambiguity,  i.e., 
with  the  rank  of  the  array  manifold  matrix.  According  to  Proukakis  and  Manikas  [27],  “any  subsequent 
research  to  handle  the  ambiguity  problem  was  mainly  concerned  with  either  the  performance  of  specific 
arrays,  or  with  the  identification  of  array  structures  free  of  ambiguities  up  to  a  certain  rank  of  ambiguity.” 

The  linear  dependence  among  the  NLA  steering  vectors  does  not  necessarily  imply  nonidentifiability,  i.e., 
even  though  the  conventional  subspace-based  DOA  estimation  algorithms  may  fail  to  deliver  unambiguous 
DOA  estimates  in  such  cases,  there  exists  other  techniques  that  could  be  developed  or  be  complemented 
MUSIC,  etc.  that  are  capable  of  manifold  ambiguity  resolution.  It  appears  that  the  failure  of  subspace- 
based  (MUSIC- type)  algorithms  to  correctly  identify  sources  with  linearly  dependent  steering  vectors  (i.e., 
to  resolve  manifold  ambiguity)  has  evolved  into  the  misunderstanding  that  “the  use  of  arrays  with  linearly 
independent  steering  vectors  is  crucial  to  applications  that  involve  estimation  of  the  DOA’s  of  multiple 
sources”  [28].  Critically  important  results  in  this  regard  were  recently  obtained  by  Proukakis  and  Manikas 
[27],  where  it  was  proven  that  practically  all  NLA  geometries  mentioned  in  the  literature  are  manifoldly 
ambiguous.  They  demonstrated  that  a  sufficient  condition  for  the  presence  of  manifold  ambiguities  in 
any  linear  array  is  that  the  antenna  array’s  aperture  measured  in  half  wavelengths.  Moreover,  in  [27]  the 
authors  introduced  an  approach  to  calculate  the  “ambiguous  generator  sets”  (AGS’s)  of  any  NLA,  i.e.,  the 
sets  of  DOA’s  that  correspond  to  linearly  dependent  steering  vectors. 

Clearly,  it  is  important  to  note  that  even  if  a  scenario  is  identifiable,  any  particular  DOA  estimation 
algorithm  may  or  may  not  yield  unambiguous  DOA  estimates.  Conversely,  if  a  particular  algorithm  (e.g., 
MUSIC)  is  applied  to  a  theoretically  identifiable  scenario,  it  may  have  the  potential  to  create  ambiguity, 
which  may  well  be  avoided  by  some  other  algorithm.  As  a  result,  one  needs  to  carefully  specify  identifiability 
and  ambiguity  aspects  with  respect  to  a  DOA  estimation  problem  and  to  clarify  and  specify  the  association 
between  manifold  ambiguity  and  nonidentifiability.  In  cases  where  manifold  ambiguity  does  not  result  in 
nonidentifiability,  one  must  attempt  to  find  an  appropriate  technique  that  can  properly  resolve  manifold 
ambiguity.  Given  that  NLA’s  are  suitable  for  DOA  estimation  beyond  the  limit  of  the  “conventional” 
number  of  independent  sources  (that  is  1  <  N  <  M),  one  needs  to  be  concerned  with  the  issue  of 
identifiability  in  the  “superior  case”  (that  is  N  >  M),  where  manifold  ambiguity  loses  its  meaning. 

In  [24],  an  effective  solution  to  the  above  problems  has  been  proposed  where  they  attempt  to  find 
a  best  fit  among  the  set  of  estimated  spatial  covariance  lags  and  source  powers  for  each  of  the  MUSIC 
DOA  estimates,  including  the  ambigious  ones.  The  proposed  approach  is  computationally  efficient  as  it 
uses  a  linear  programming  algorithm  and  has  demonstrated  a  high  probability  of  correct  identification 
in  manifoldly  ambiguous  scenarios  exceeding  the  performance  of  “abnormal”  MUSIC  behavior  that  is 
well-known  and  caused  by  the  presence  of  low  SNR  and/or  number  of  data  samples. 

One  of  the  NLA  subclasses  considered  is  the  “fully  augment  able”  arrays,  which  have  the  property  that 
all  intermediate  distances  (lags)  are  available,  that  is  there  is  no  gaps  in  the  set  of  inter-sensor 

separations  implying  that  the  set  of  all  distances  is  complete.  The  problem  of  DOA  estimation  for  the 
fully  augmentable  sparse  linear  arrays  when  there  are  more  uncorrelated  sources  than  sensors  in  a  NLA 
has  been  investigated  in  [29]. 

The  above-mentioned  algorithms  apply  almost  exclusively  to  linear  arrays;  whereas  in  most  of  the 
real-life  situations  the  physics  of  the  problem  dictates  a  non-linear,  e.g.  circular  or  even  volumetric,  array 
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geometry.  Additionally,  these  cannot  be  used  in  cases  where  the  array  is  sparse  leading  to  spatially  aliased 
data.  At  present,  there  are  only  a  few  DOA  estimation  approaches  for  volume  arrays  (Liu  &  Mendel;  1998 
[30]),  (Barthelemy  &  Willett;  1998  [31])  and  for  non-uniform  spatial  sampling  (Zoltowski  &  Mathews;  1994 
[32])  but  their  applicability  to  multiple  target  detection  and  tracking  is  doubtful.  This  is  particularly  true 
when  the  targets  are  spatially  close  together  and  in  presence  of  heavy  clutter  and  correlated  noise. 


In  general,  for  a  non-uniform  volume  array,  the  time  delays  can  be  related  to  the  DOA  of  the  source 
and  the  actual  location  of  the  sensors  i.e. 
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where  Kn  =  [ sinOncosipn ,  sinOnsimpn ,  cos0n]T  is  the  DOA  vector  for  the  nth  source  with  8n  and  y?n  being 
the  azimuth  and  elevation  angles,  respectively,  Xm  =  [ximiX2m,X3rn]T  represents  the  coordinate  vector  of 
the  mth  sensor.  In  this  case,  the  steering  matrix  A  becomes, 


A(9,<p)  =  [a(6 i,<pi),-  ■  * 

with  the  steering  vector  a(0n,  (pn)  =  [1  e~iUJcT2n,  •  •  •  ,  e~JWcTMn]T. 


Thus,  in  addition  to  the  difficulties  in  resolving  multiple  closely  spaced  targets  another  challenge  in  this 
development  is  how  to  detect  and  extract  useful  target  information  in  spatially  aliased  data  due  to  sparse 
nature  of  most  of  the  sensor  arrays  used  for  this  problem.  This  is  especially  important  in  presence  of  array 
manifold  ambiguity. 


2.4  DOA  Estimation  for  Multiple  Invariance  Arrays 

The  main  limitation  for  the  MI-ESPRIT  is  that  it  requires  a  multidimensional  search  for  the  DOA  arrays, 
which  is  computationally  very  time-consuming.  In  (Swindlehurst,  Stoica,  &  Jansson;  2001  [4])  a  new 
algorithm  is  proposed  that  attempts  to  generalize  the  MUSIC  technique  denoted  as  MI-MUSIC.  Unlike 
MI-ESPRIT  the  DOA  estimate  is  obtained  only  through  a  1-D  search  as  opposed  to  a  multi-dimensional 
search,  however  unlike  ESPRIT  it  does  not  enjoy  the  statistical  optimality  properties.  Furthermore,  the 
authors  have  also  generalized  the  Root-MUSIC  and  MODE  techniques  to  multiple  invariance  versions. 
The  advantage  of  the  above  algorithms  is  that  even  the  need  for  a  1-D  search  is  eliminated.  It  has  been 
shown  that  MI-MUSIC  and  MI-Root-MUSIC  algorithms  yield  good  performance  and  that  MI-MUSIC  may 
actually  be  more  robust  than  MI-ESPRIT  in  difficult  threshold  scenarios.  For  highly  correlated  sources 
both  MI-ESPRIT  and  MI-MODE  perform  better  than  MI-MUSIC  and  MI-Root-MUSIC.  The  details  are 
briefly  stated  below. 


In  the  multiple  invariance  array  [4]  it  is  assumed  that  we  have  p  identical  translated  subarrays  of  m 
elements  each  that  are  arranged  with  arbitrary  but  known  spacing.  In  order  to  describe  the  array  formally 
one  of  the  subarrays  is  considered  as  the  reference  and  the  known  translational  and  angular  displacements 
of  the  other  subarrays  are  denoted  by  8k  and  ak,  for  k  =  1,  •  •  •  ,p- 1.  The  structure  may  now  be  expressed 
as 


JA{6)  =  A 


F 


p- 1 


where  J  is  a  known  mp  x  M  selection  matrix,  F  is  the  unknown  response  of  the  reference  subarray,  and 
is  a  diagonal  matrix  defined  by 

=  diag{ej*sm{0l+ak),  ■■■  ,  eP™n(ei+olk)} 
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where  8k  is  measured  in  the  half-wavelengths  units.  If  the  subarrays  are  non-overlapping  then  the  rows  of 
A(9)  may  be  arranged  so  that  J  =  I.  Although  the  condition  d  <  M  is  necessary  for  any  DOA  estimation 
algorithm  based  on  the  covariance  matrix,  it  turns  out  that  it  may  not  be  sufficient  in  our  application 
given  that  there  are  additional  parameters  that  need  to  be  estimated  in  F.  Briefly  stated  the  MI-ESPRIT 
algorithm  applied  to  our  configuration  amounts  to  the  following  procedure 

F,0,T  =  arg  minF,$,T\\JGsW2  -  AT\\2F 

where  ||.|| F  denotes  the  Frobenius  norm,  and  IF  is  a  suitably  chosen  weighting  matrix,  and  T  is  a  full-rank 
matrix  of  dimension  d  x  d  that  satisfies  JGS  =  AT.  The  necessary  condition  for  the  multiple  invariance 
array  is  that  d  <  min{M  —  l,rn(p  —  1)}.  Similarly,  the  extensions  for  the  algorithm  to  MI-MUSIC  and 
MI- Root-MUSIC  may  also  be  developed  as  in  [33]. 

2.5  DOA  Estimation  for  Non-Gaussian  Sources 

Various  DOA  estimation  algorithms  are  developed  based  on  the  covariance  matrix  of  the  array  output. 
The  delay-and-sum  method  as  well  as  Capon’s  minimum  variance  method  use  the  sample  covariance 
matrix  to  measure  the  power  as  a  function  of  the  DOA.  The  subspace  methods  (such  as  MUSIC,  ESPRIT, 
minimum  norm,  etc.)  estimate  the  signal  and  noise  subspaces  from  the  eigenvectors  of  the  covariance 
matrix.  Provided  that  the  data  is  Gaussian  distributed,  the  sample  covariance  matrix  is  an  optimal 
estimator  of  the  unknown  covariance  matrix.  However,  as  it  is  the  case  in  most  real-life  situations  the 
data  are  generally  non-Gaussian,  which  can  result  in  poor  performance  and  unreliable  DOA  estimate  of 
these  algorithms.  In  another  class  of  the  DOA  algorithms  that  are  based  on  the  maximum  likelihood  (ML) 
principle,  the  estimates  for  the  DOA’s  are  derived  by  assuming  that  the  noise  is  modeled  as  a  Gaussian 
random  process.  Similar  to  the  previous  algorithms,  under  these  circumstances  the  ML  estimates  are 
significantly  under-performing  and  are  highly  unreliable.  In  other  words,  when  the  sensor  noise  is  spatially 
non-uniform  white  process,  neither  the  conventional  ML  methods  nor  the  colored  noise  modeling  ML 
techniques  are  expected  to  provide  acceptable  performance. 

In  the  general  case,  the  sensor  noise  must  be  treated  as  a  colored,  spatially  dependent  process.  For 
example,  the  ambient  noise  and  wind  conditions  have  been  observed  to  be  clearly  non-Gaussian.  This  obvi¬ 
ously  calls  for  and  necessitates  development  of  “robust”  DOA  estimation  algorithms  for  these  nonstandard 
circumstances.  Current  research  in  the  literature  attempting  to  address  the  above  issues  are  reported  in 
(Lee  &  Kashyap;  1992  [34])  using  M-e stimation.  The  problem  of  DOA  estimation  in  heavy-tailed  noise 
generated  from  an  ct-stable  distribution  have  also  been  studied  in  (Tsakalides  &  Nikias;  1996  [35]).  Some 
authors  attempt  to  model  the  noise  as  a  finite  mixture  of  Gaussian  random  variables  and  use  the  SAGE 
algorithm  as  in  (Kozick  &  Sadler;  2000  [36]). 

The  above-mentioned  robust  algorithms  for  the  DOA  estimation  all  suffer  from  the  limitations  that  they 
rely  on  the  knowledge  of  the  pdf  or  the  number  of  the  mixtures.  Furthermore,  they  also  need  user-specified 
parameters  (such  as  thresholds  and  weight  matrices)  that  have  to  be  selected  based  on  trial- and-err or  or 
heuristic  methods,  lending  these  algorithms  somewhat  sensitive  to  the  initialization  problems.  Therefore, 
there  is  a  strong  interest  in  developing  computationally  practical  and  easy  to  implement  DOA  estimation 
algorithms  in  real-time  for  the  non-uniform  noise  situation.  In  (Pesavento  &  Gershman;  2001  [37])  this 
problem  is  tackled  by  developing  a  deterministic  ML  DOA  estimator  based  on  a  non-uniform  noise  model. 
The  algorithm  is  implemented  based  on  the  so-called  stepwise  concentration  of  the  log-likelihood  function 
with  respect  to  the  signal  and  noise  nuisance  parameters  and  requires  only  a  few  iterations  to  converge. 
In  (Visuri,  Oja  and  Koivunen;  2001  [38])  a  high-resolution  DOA  estimation  algorithm  is  proposed  that  is 
based  on  estimating  the  signal  and  noise  subspaces  from  spatial  sign  or  rank  covariance  matrix.  Convergent 
estimates  of  the  non-parametric  statistics  are  developed.  Furthermore,  they  have  shown  that  their  proposed 
algorithms  have  comparable  performance  compared  to  the  optimal  methods  in  presence  of  Gaussian  noise 
and  are  not  dependent  on  user-specified  parameters. 
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2.6  DOA  Estimation  Under  Missing  Data/ Array  Scenarios 

The  problem  of  estimating  the  DOA  from  incomplete  or  partially  available  measurements  has  recently 
been  investigated  by  (Larsson  &  Stoica;  2001  [39]).  There  are  many  factors  that  could  cause  faults  in 
arrays  which  could  eventually  lead  to  failed  arrays,  resulting  in  incomplete  available  information.  In  our 
particular  problem  such  failure  may  be  resulted  from  the  environmental  and  natural  conditions.  The 
presently  available  DOA  techniques  are  not  designed  to  cope  and  deal  with  these  uncertainties,  i.e.,  they 
are  not  “fault-tolerant”.  Specifically,  the  techniques  that  use  the  sample  covariance  matrix  as  an  estimate 
of  the  array  output  covariance  matrix  could  fail  to  yield  accurate  and  reliable  results  as  they  are  based  on 
the  assumption  that  the  sample  covariance  matrix  is  an  accurate  representation  of  the  array  covariance 
matrix. 

The  model  of  the  sensor  array  that  was  considered  in  [36]  assumes  that  one  or  more  sensors  fail  to 
work  after  a  certain  period  of  time.  Specifically,  assuming  that  the  total  number  of  sensors  in  the  array 
is  Mi  and  that  data  were  collected  with  all  M\  sensors  for  t  =  1,  •  •  •  ,  L\.  After  t  =  Li,  at  least  one 
sensor  fails  to  work  so  that  M2  working  sensors  still  remain.  Using  these  sensors  data  is  collected  from 
t  —  Li  +  1,---  ,  L\  +  L2-  This  process  is  repeated  for  t  >  L\  +  L2.  It  is  clear  that  Mi  >  M2  >  •  *  >  Mp. 
It  is  assumed  that  {Mi,  •  •  •  ,  Mp}  and  {Li,  •  ■  •  ,  Lp}  are  known  a  priori .  This  is  justified  by  the  fact  that 
failing  sensors  and  the  times  when  they  stop  functioning  are  typically  easy  to  identify. 

The  contribution  of  the  work  by  (Larsson  &  Stoica;  [39])  is  to  enable  the  existing  techniques  such 
as  MUSIC,  ESPRIT,  etc.  with  a  reliable  statistical  estimate  of  the  array  covariance  matrix  to  be  used 
instead  of  the  sample  covariance  matrix.  Several  techniques  were  considered.  In  the  first  approach,  the 
ML  estimate  of  the  covariance  matrix  and  its  asymptotic  accuracy  were  derived  and  discussed.  It  was 
shown  that  by  using  the  Covariance  Matching  technique  (COMET),  as  opposed  to  MUSIC,  one  can  use 
information  on  the  accuracy  of  the  covariance  matrix  estimate,  and  that  this  information  can  improve 
the  performance  in  estimating  the  DOA’s  from  the  incomplete  data.  It  was  concluded  that  the  proposed 
COMET  approach  may  be  the  preferred  choice  when  treating  parameter  estimation  from  incomplete  data. 

2.7  DOA  Estimation  using  Independent  Component  Analysis  and  Blind  Source  Sep¬ 
aration 

Blind  source  separation  (BSS)  is  an  approach  used  to  estimate  original  source  signals  through  only  the 
information  of  the  mixed  signals  observed  in  each  input  channel.  Recently  work  has  been  conducted  for 
the  BSS  based  on  the  independent  component  analysis  (ICA)  [40].  There  are  several  methods  proposed  in 
the  literature  in  which  the  inverse  of  the  complex  mixing  matrices  are  calculated  in  the  frequency  domain 
to  deal  with  the  arrival  lags  among  each  of  the  elements  of  the  microphone  array  system  [41,42].  However, 
this  ICA-based  approach  has  the  disadvantage  that  there  is  difficulty  with  the  low  convergence  of  nonlinear 
optimization  [43].  In  [44]  a  new  algorithm  is  proposed  for  BSS,  in  which  ICA  and  beamforming  are  combined 
to  resolve  the  low-convergence  problem  through  optimization  in  ICA.  The  proposed  method  consists  of  the 
following  three  parts:  (1)  frequency-domain  ICA  with  DOA  estimation,  (2)  null  beamforming  based  on  the 
estimated  DOA,  and  (3)  integration  of  (1)  and  (2)  based  on  the  algorithm  diversity  in  both  iteration  and 
frequency  domain.  The  inverse  of  the  mixing  matrix  obtained  by  ICA  is  temporally  substituted  by  the 
matrix  based  on  null  beamforming  through  iterative  optimization,  and  the  temporal  alternation  between 
ICA  and  beamforming  can  realize  fast-  and  high-convergence  optimization.  The  results  of  the  signal 
separation  experiments  have  shown  that  the  signal  separation  performance  of  the  proposed  algorithm  is 
superior  to  that  of  the  conventional  ICA-based  BSS  method,  even  under  reverberant  conditions. 

2.8  Multiple  Target  Tracking  and  DOA  Estimation  for  Unattended  Ground  Sensors 

The  work  at  the  Army  Research  Laboratory  (ARL)  (Pham  &  Fong;  1997  [2]),  (Pham  &  Sadler;  1997 
[3])  involved  development  of  array  signal  processing  algorithms  using  small  baseline  acoustic  arrays  to 
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estimate  DOA’s,  track,  and  classify  moving  targets  from  their  acoustic  signatures.  Since  conventional 
beamforming  schemes  generally  do  not  provide  good  results  in  these  situations,  owing  to  the  structural 
limitations  of  the  array  and  lack  of  spatial  coherence,  high-resolution  DOA  algorithms  must  be  developed. 
Current  efforts  on  incoherent  and  coherent  wideband  MUSIC  show  some  promise  for  detecting  and  tracking 
multiple  ground  vehicles.  In  [2]  by  exploiting  the  multi-spectral  content  of  the  targets,  better  results  were 
reported  using  wideband  signal  subspace  DOA  algorithms.  Experimental  results  using  a  circular  array 
of  six  sensors,  with  a  diameter  of  8  ft  indicated  the  advantages  of  the  wideband  MUSIC  over  the  delay- 
sum  (DS)  algorithm.  In  [3]  focused  wideband  adaptive  array  processing  algorithms  are  employed  for  the 
high-resolution  DOA  estimation  of  ground  vehicles.  Several  methods  including  steered  covariance  matrix 
(STCM)  [45]  and  spatial  smoothing  or  interpolation  using  the  array  manifold  interpolation  (AMI)  are 
considered.  In  addition,  this  study  provides  experimental  analysis  of  incoherent  and  coherent  wideband 
MUSIC  algorithms  and  narrowband  MUSIC  for  a  circular  array  of  acoustic  sensors  in  terms  of  accuracy  and 
stability  of  DOA  estimates.  It  was  shown  that  for  certain  SNR  conditions  incoherent  wideband  methods 
yield  more  accurate  DOA  estimates  than  the  coherent  methods  for  highly  peaked  spectra  while  for  sources 
with  flat  spectra  the  coherent  wideband  methods  generate  more  accurate  DOA  estimates.  The  study  are 
indicated  that  the  coherent  wideband  methods  are  much  more  statistically  stable  than  the  incoherent  ones. 
The  price  paid  for  this  improved  stability  is  high  computational  cost. 

More  recently,  (Boettcher,  Sherman  &  Shaw;  2002,  [46])  developed  a  simple  time-difference  of  arrival 
(TDOA)  estimation  method  for  localization  of  targets  using  a  network  of  acoustic  sensor  arrays.  TDOA  is 
performed  using  spectral  peaks  from  the  acoustic  signatures  to  allow  drastic  reduction  in  the  bandwidth 
required  to  collaboratively  determine  bearing  angles.  Although  this  scheme  shares  some  of  the  advantages 
of  the  coherent  and  closest  point  of  approach  (CPA)-based  algorithms,  it  cannot  be  used  to  estimate 
DOA’s  associated  with  multiple  closely  spaced  targets  or  in  difficult  target  scenarios.  In  [47],  Kaplan,  Le 
and  Molnar  (2001)used  ML  method  to  localize  a  moving  target  using  a  network  of  passive  acoustic  arrays 
where  each  node  consists  of  an  array  of  microphones  configured  in  a  “wagon  wheel”  pattern.  Each  node 
transmits  a  DOA  estimate  to  a  central  processor.  The  central  node  processes  the  bearing  information  to 
localize  and  track  targets  of  interest.  A  distributed  sensor  measurement  strategy  dynamically  designates 
the  active  nodes  and  assigns  nodes  to  central  processing  task. 


2.9  Target  Classification  Based  Upon  Acoustic  Signatures 

Although  there  are  a  number  of  bearing  angle  estimation  methods  for  detection  and  tracking  of  vehicles 
(wheeled  and/or  tracked)  using  UGS  acoustic  arrays,  the  literature  on  classification  of  these  targets  based 
upon  their  acoustic  signature  and  other  features  is  very  limited.  Wu,  Siegel  and  Khosla  [48]  used  Principal 
Component  Analysis  (PCA)  of  the  windowed  acoustic  signatures  in  the  frequency  domain  as  features  for 
classification.  When  a  new  windowed  signal  is  applied,  its  spectrum  vector  is  projected  onto  the  principal 
eigenvectors.  Results  on  discrimination  of  passenger  cars  against  trucks  are  also  presented. 

We  believe  subband  features  as  well  as  velocity,  acceleration  and  curvature  profiles  of  various  detected 
targets  can  be  used  to  discriminate  between  different  types  of  targets.  The  subband  features  are  expected 
to  reveal  subtle  tonal  features  in  different  frequency  bands  and  provide  clues  to  classify  different  heavy  and 
light  vehicles  of  both  wheeled  and  tracked  types.  Of  particular  importance  is  to  determine  how  the  subband 
features  can  aid  in  extracting  information  about  light  vehicles  as  they  present  a  major  challenge  to  both 
detection  and  classification  systems.  The  developed  classification  system  must  be  flexible  in  that  other 
target  clues  from  other  types  of  sensors,  e.g.  seismic,  which  complements  the  decision  making  process,  may 
easily  be  incorporated  into  the  system.  The  subband  acoustic  features  together  with  the  targets  dynamical 
features  can  be  measured  at  every  few  sensor  snapshots  and  then  used  as  inputs  to  a  target  classification 
system.  The  sequential  Bayes  classification  scheme  [49]  developed  by  the  PI  can  be  employed  here.  The 
results  can  be  updated  on-line  as  more  data  is  becoming  available. 
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3  Generating  Synthesized  Multiple  Target  Scenarios 

Two  data  sets  were  received  from  the  US  Army  TACOM-ARDEC  in  Picatinny  Arsenal,  NJ.  The  first  data 
set  (SAFE  II)  was  collected  using  three  nodes  each  consisting  of  a  circular  (wagon  wheel-type  pattern)  array 
of  five  identical  elements.  The  second  data  set  (YUMA)  was  collected  using  four  nodes  with  a  similar  array 
structure.  Although  these  data  sets  contain  multiple  target  scenarios,  according  to  the  technical  point  of 
contact  (POC),  Mr.  Robert  Wade,  they  are  not  representative  of  real  battlefield  scenarios  where  multiple 
closely  spaced  targets  in  abreast,  staggered  or  single-file  formation  are  encountered.  Consequently,  several 
methods  for  generating  synthetic  multiple  target  data  sets  were  developed  and  investigated.  These  data 
sets  are  then  used  to  develop  and  test  high-resolution  and  high  confidence  DOA  estimation  and  source 
separation  algorithms  that  offer  reduced  confusion  in  separating  multiple  targets  in  difficult  battlefield 
situations.  The  following  subsections  present  two  different  methods  for  generating  synthesized  data  sets. 

3.1  Autoregressive  (AR)-Based  Method 

In  these  simulations  the  number  of  snapshots  is  set  to  N  =  1024  and  the  sampling  frequency  fs  =  1024  Hz. 
The  number  of  sensors  is  L  =  5  and  the  number  of  sources  is  P  >  2  unless  otherwise  specified.  Two  types 
of  sources  are  used.  The  sources  are  generated  using  1st  order  Autoregressive  (AR)  processes  driven  by 
mutually  independent  zero-mean  white  Gaussian  noise  with  unit  variances.  That  is  for  the  ith  source  we 
can  write 


Si(n)  =  a,iSi(n  -  1)  4-  Vi(n ) 


(1) 


where  Vi(n )  is  a  zero-mean  white  Gaussian  noise  with  unit  variance  mutually  independent  of  Vj(n),j  7^  i. 
Variance  of  each  source  is  adjusted  to  unit  before  the  formation  of  the  array  measurements.  Some  of  the 
truth  files  associated  with  SAFE  II  data  were  used  to  generate  the  synthetic  target  tracks  with  various 
close  formations.  Ten  (10)  frequency  bins  ( h  =  10)  are  used  in  the  wideband  DOA  estimation  algorithms 
(As  described  in  details  in  Section  4).  Bandpass  filters  had  to  be  used  to  achieve  desired  frequency  bands 
for  the  sensor  measurement  in  the  STCM  and  WSF  methods  (refer  to  Section  4  for  details).  These  are  FIR 
type  filters  with  21  weights  designed  using  the  appropriate  MATLAB  function.  The  frequency  range  was 
selected  within  the  interval  [/l,/h]  =  [200  Hz,  300  Hz],  Simulations  with  and  without  the  passband  filter 
were  conducted  for  16  independent  runs.  In  the  case  of  the  circular  array  geometry  the  sensor  spacing  dx 
and  dy  are  set  to  be  2ft.  («  0.6  (m))  as  in  the  SAFE  II  data. 

The  number  and  selection  of  the  frequency  bins  influence  the  performance  of  the  algorithms  to  some 
extend,  and  its  choice  needs  to  be  decided  carefully.  The  selection  of  the  frequency  bins  is  an  issue  that 
would  need  further  research  to  optimize  the  performance  of  the  algorithms.  In  our  simulations,  the  powers 
of  all  the  sensors  are  put  together  at  each  frequency  bin,  and  the  frequency  bins  that  indicate  the  largest 
values  are  selected  one  by  one  until  the  desired  number  is  reached.  If  the  bin  number  is  too  small,  the 
frequency  bins  that  contribute  to  the  accurate  DOA  estimate  may  be  missed.  On  the  other  hand,  if  the 
bin  number  is  too  large,  the  heavily  noisy  frequency  bins  may  be  included,  leading  to  poor  DOA  estimates. 
This  ad-hoc  way  of  selection  works  rather  well  in  the  simulations,  and  the  number  h  =  10  was  found  to  be 
a  good  choice  in  these  results. 

3.2  Interpolation-Based  Method 

More  realistic  simulated  data  sets  containing  multiple  targets  in  abreast,  single-file  and  staggered  formations 
can  be  generated  using  an  interpolation-based  method  that  accounts  for  Doppler  effects  off  moving  sources. 
This  method  generates  substantially  more  accurate  simulated  data  sets  as  it  does  not  assume  time  delays 
are  multiples  of  the  sampling  interval.  The  wagon-wheel  circular  array  geometry  with  4ft.  radius  was 
considered  and  tested  on  the  generated  data  sets. 
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To  understand  the  difference  between  the  previous  approach  and  this  new  method  imagine  two  hypo¬ 
thetical  scenarios.  In  first  case,  assume  that  the  track  that  the  vehicle  traverses  is  known  and  a  number  of 
signal  emitters  are  placed  on  the  track  at  a  fixed  interval.  These  emitters  emit  the  source  signal  and  are 
turned  on/off  by  a  remote  controller.  Using  the  original  position  of  the  virtual  target,  and  its  velocity,  one 
can  compute  the  time  instances  at  which  the  target  passes  those  positions.  Then,  to  synthesize  the  data, 
one  can  use  the  remote  controller  to  turn  on  the  emitters  one  by  one  at  these  instances  and  collect  the 
source  signals  using  the  array  of  sensors/microphones.  Most  of  the  commonly  used  schemes,  including  the 
above-  mentioned  AR-based  method,  for  generating  simulated  array  data  are  based  upon  such  simplified 
assumption.  In  the  second  hypothetical  scenario,  assume  that  instead  of  placing  a  large  number  of  emitters 
along  the  track,  only  one  signal  emitter  is  placed  on  an  ideally  silent  vehicle  which  will  follow  the  track 
with  the  same  velocity  as  the  target.  Clearly,  comparing  these  two  scenarios,  one  can  see  that  the  latter 
scheme  has  several  advantages.  The  first  scheme  can  only  produce  the  data  at  some  isolated  time  intervals, 
while  using  the  second  scheme  accurate  source  signals  are  generated.  Most  importantly,  the  first  scheme 
does  not  take  into  account  for  the  Doppler  effects,  whereas  in  second  scheme,  this  is  built  in  the  setup.  In 
practice,  what  we  can  assume  to  know  is  the  signal  emitted  from  a  stationary  target  that  is  sampled  at 
certain  sampling  frequency. 

In  our  simulations,  we  assume  that  the  signals  from  different  targets  and  noise  are  additive  at  the 
sensors.  Under  this  assumption,  we  can  first  calculate  the  signals  received  by  a  sensor  from  different 
targets  separately  in  a  one-target  one-sensor  fashion.  Then,  all  the  signals  from  different  targets  are  added 
together  with  some  noise.  Figure  1  illustrates  this  one-target-one-sensor  scenario.  The  position  of  the 
moving  target  at  time  t0  is  denoted  by  (x(t0),  y{to))  and  x(t ),  y(t)  at  time  t .  The  position  of  the  stationary 
sensor  is  ( xs,ys ).  The  velocity  of  the  target  is  vy(t)),  where  vx(t)  and  vy(t)  denote  the  velocity 

components  along  the  x-axis  and  y-axis,  respectively.  The  velocity  of  the  target  could  be  a  time  varying 
function. 


v,<t)  v(t) 


Figure  1:  The  scenario  of  one-target-one-sensor. 
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(2) 


At  time  t,  t  >  to?  the  position  of  the  target  is  given  by 

x(t)  —  x(to)  +  f  vx(t)dt 

Jt0 

y(t)  =  y(t0)+  [  vy(t)dt  (3) 

Jto 

Thus,  the  distance  between  the  target  and  the  sensor  is  d{t)  =  \/{x(t)  —  xs\2  +  [y2(t)  —  ys}2.  Let  us 
assume  that  signal  5(f)  emitted  by  the  target  is  known  in  advance  and  let  c  be  the  speed  of  sound  travelling 
in  air.  Then,  the  signal  received  by  the  sensor  is  R(t)  =  5(f  -  d(t)/c).  Note  that  this  model  actually  takes 
into  account  the  Doppler  effects  caused  by  the  movement  of  the  target.  This  will  be  more  evident  if  we 
assume  that  the  target  is  moving  directly  toward  the  sensor  with  a  constant  speed,  the  velocity  in  this 
case  can  be  written  as  (—  x(t2^xdf  _  wjiere  j'  js  an  arbitrary  interval.  Then  the  resultant  signal 

received  by  the  sensor  is  R(t)  =  5((1  +  ^l)f  —  It  is  obvious  that  the  term  1  +  multiplied  by  t 

indicates  that  R(t)  is  a  frequency  shifted  version  of  5(f) . 

In  practice,  we  assume  that  the  sampled  version  of  the  wideband  signal  5(f)  emitted  by  the  target,  i.e. 
s(n)  =  S(nTo )  where  To  is  the  sampling  interval,  can  be  modelled  by  an  autoregressive  (AR)  model.  It 
is  essential  that  when  we  simulate  the  signal  r(n)  received  by  the  sensor,  we  also  use  the  same  sampling 
interval  Tq.  Then  we  have 


r(n)  =  R(nTo) 

—  S(nTo  -  d(nTo)/c )  —  S(tr) 


where  tr  :=  tiTq  —  d(nTo)/c.  The  problem  is  that  tT  is  not  always  a  multiple  of  To,  which  precludes  the 
direct  use  of  values  of  s(n).  The  solution  is  to  use  interpolation.  In  this  case,  let  ty  and  tc  be  the  floor  and 
ceiling  of  tr,  respectively,  i.e.  f/  =  |_Aj  and  tc  =  [tr] .  Note  that  tc  —  tf  =  To,  then 


r(n) 


To 


(tr  tf )  +  S(tf ) 


To  compute  d(nT0),  both  x(tiTq)  and  y(nTo)  are  needed.  Equations  (2), (3)  can  be  approximated  as 
follows: 

n 

x(n%)  =  X'(f0)  +  y^vx{{i  -  1)T0)T0 

i= l 
n 

y(nT0)  =  y(t0 )  +  -  1)T0)T0 

2=1 

Figure  2  illustrates  the  positions  of  the  array  nodes  and  targets  in  the  case  where  two  abreast  (side-by- 
side)  targets  are  present.  The  two  targets  are  moving  with  constant  velocities  and  hence  remain  abreast  in 
all  observations.  The  signals  received  by  the  sensors  are  generated  separately  and  packed  into  a  database 
for  later  usage. 

To  generate  realistic  source  signals,  we  devised  a  scheme  where  a  5th  order  AR  model  is  fitted  to  the 
source  signal  in  one  of  the  real  SAFE  II  data  files  that  contained  only  one  target.  The  AR  model  is  found 
to  be 


x(n)  -  1.627x(n  -  1)  -  1.117a;(n2)  +  0.6719rr(n  -  3)  -  0.4124x(n  -  4)  +  0.1698x(n  -  5)  +  7?(n) 


12 


Figure  2:  The  setup  of  the  case  of  two  abreast  targets. 

where  77(71)  is  a  zero- mean  white  Gaussian  process  with  variance  a2n  =  145.28. 

Figure  3(a)  illustrates  an  example  of  two  simulated  target  tracks  generated  using  this  scheme  where  the 
two  targets  are  far  apart.  Figure  3(b)  shows  the  detected  tracks  using  the  STCM  method.  This  method 
is  described  in  details  in  Section  4.1.  Only  one  array  is  used  in  this  DO  A  estimation.  We  can  see  that  the 
two  targets  are  perfectly  detected,  which  suggests  that  the  synthesized  data  is  indeed  useful. 

4  High-Resolution  Wideband  DOA  Estimation  Methods 

One  of  the  main  issues  of  paramount  importance  in  this  Phase  I  research  project  was  the  development 
of  a  high-resolution  wideband  DOA  estimation  scheme  for  separating  multiple  closely  spaced  targets.  We 
investigated  the  applicability  of  several  DOA  estimation  methods  for  the  problem  in  hand.  Three  different 
wideband  DOA  estimation  methods,  namely  incoherent  weighted  MUSIC  (IWM)  [33],  coherent  steered 
covariance  matrix  (STCM)  [2], [45], [50]  and  weighted  subspace  fitting  (WSF)  [51]  are  considered.  These 
algorithms  are  known  in  array  signal  processing,  but  many  aspects  of  their  performance  are  not  completely 
understood  and  investigated,  especially  when  applied  to  unattended  passive  acoustic  arrays.  First,  these 
three  algorithms  are  briefly  introduced  in  a  way  that  facilitates  their  implementation.  More  theoretical 
details  can  be  found  in  the  related  references.  Simulations  on  both  real  (SAFE  II)  and  synthesized  data  are 
conducted  to  compare  the  performance  of  these  algorithms,  their  strengths  and  weaknesses  and  to  identify 
the  issues  that  need  to  be  carefully  considered  when  these  are  applied  to  our  problem.  Extensive  simulation 
results  are  presented  for  the  uniform  linear  array  (ULA)  as  well  as  the  wagon-wheel  array.  Finally,  some 
key  observations  are  provided. 
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Track  of  real  target  detected  DOA 


Figure  3:  The  real  and  detected  tracks  of  two  far- apart  targets. 


4.1  Steered  Covariance  Matrix  (STCM)  Algorithm 

In  this  section  we  briefly  review  the  STCM  method  and  present  some  preliminary  results  on  both  the 
synthesized  and  real  SAFE  II  data. 


Let  us  denote  the  complex  analytic  signal  of  the  output  of  a  sensor  at  location  xm  by  y(ra,£),  for 
m  =  —  K, . . . ,  K)  where  the  total  number  of  sensors  is  M  =  2 K  +  1.  Then,  y(m,  t ),  observed  over  a  time 
interval  of  T  seconds  can  be  represented  as, 

P 

y(m,t )  =  ^2si[t-  Tm(0i)]  +  v(m,t) 

i—1 

where  Si(t)yi  =  1, . . . ,  P  are  stationary,  zero-mean  processes  corresponding  to  the  received  source  signals 
and  is  the  spatially  and  temporally  stationary  zero-mean  noise  process  at  the  mth  sensor.  The 

second-order  statistics  of  the  field  are  completely  specified  by  the  cross-spectral  density  matrices,  R(u>k)  = 
E[Y(ujk)Y(u)k)*],foi  k  —  Z,..,h,  where  E[]  denotes  the  expectation  operator  and  Y(cvk)  is  the  frequency 
component  at  frequency  cjk  .  For  the  above  model,  we  have 


R(ojk)  =  A(vk,6)Ps(u)k)A(u)k,0y  +  Rv(ujk) 


where  A(ujkj  0)  =  [a(cuk,  Qi),  a(uk,  62), . . . ,  a(u>ky  0p)]  is  the  M  x  P  source  direction  matrix,  a(ujky6i)  ~ 
[e-jukT-K(ei) ^  ^  Js  direction  vector  of  the  ith  source,  Ps( u>k)  is  a  P  x  P  source  spectral 

density  matrix,  9  =  [0i, . . . ,  0p],  and  Rv(tok)  represents  the  M  x  M  noise  covariance  matrix.  For  wideband 
sources  the  problem  is  that  the  same  source  contributes  a  different  rank-one  component  to  Rv(u)k)  at  the 
temporal  frequency.  To  avoid  this  ambiguity  problem,  the  basic  idea  behind  the  STCM  [45], [50]  is  to 
preprocess  the  sensor  outputs  to  ensure  that  arrivals  from  a  particular  direction  have  the  same  rank-one 
representation  at  all  temporal  frequencies.  Letting  8  to  denote  the  direction  of  interest,  we  define  Ts(8,uk) 
as 

'  e3^kT-K{0)  Q  ...  0 

0  ejuJkT_K+l(0)  ...  0 


TsiO.ujk)  = 


(4) 


0 


0  ... 
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Now,  by  multiplying  T^O^k)  with  Y(ujk),  we  get 

h 

ys(t,8)  =  Y,Ts(6,uk)Y(uky^t  (5) 

k=l 

Substituting  (4)  into  (5)  we  get 

ys(t,9)  =  {y[-K,t  +  T-K{0)\,y{-K  +  l,t  +  T_/r+i(0)],  ■  •  • , y[K,t  +  Tk{6)]}T  (6) 

The  steered  covariance  matrix  corresponding  to  direction  6  is  given  by 

h 

Rs(8)  =  '52Ta{u>k,9)R(u>k)T,(u>k,0y  (7) 

k=l 

The  structure  of  Rs(0)  can  be  used  to  motivate  a  variety  of  wideband  spatial  spectral  estimation 
methods.  For  the  circular  array  in  our  problem,  the  STCM  is  implemented  as  follows. 

1.  Form  the  estimated  cross-spectral  density  matrices,  R(u*k)i  over  the  frequency  band  of  interest. 

A  common  method  of  forming  R(u) )  from  discrete-time  sensor  outputs  is  to  divide  the  T  second 
observation  into  N  non-overlapping  segments  of  AT  seconds  each,  and  then  apply  DFT  to  obtain 
uncorrelated  frequency  domain  vectors,  Yn(u^),  for  each  segment  n  =  1, . . . ,  N.  The  cross-spectral 
density  matrix  at  frequency  is  then  estimated  by  taking 

1  N 

R(uk)  =  -J^Yn(uJk)Yn(uky 

71=1 

2.  For  each  steering  direction,  0,  of  interest,  compute  the  estimated  steered  covariance  matrices, 

h 

Rs(6)  =  ^2Ts(cjk,9)R(ujk)Ts(cok,8)* 

fc=Z 

where  are  the  frequency  band  of  interest,  and  Ts(u>k,  0)  for  a  circular  array  of  five  elements  is 

Ts(uk,  0)  =  Diag{e~^UJkdyS:irie/c^  e-i^^cos^/c}  i;  ejukdysmO/c,  ejujkdx  cos 0/cj 

3.  Compute  i?s(0)-1  and  form  the  “Steered  Minimum  Variance”  spectral  estimator  Zstrnv{0)  —  [1*J?S(0)_11] 
where  1  is  an  M  x  1  vector  of  all  ones,  for  each  steering  direction  8  to  obtain  a  broad-band  spatial 
power  spectral  estimate.  The  estimation  of  the  source  locations  is  achieved  by  determining  the  peak 
positions  of  the  spatial  power  spectral  estimate  ZstmV(0)- 

4.2  Incoherent  Wideband  MUSIC  (IWM)  Algorithm 

In  the  array  considered  here,  we  assume  there  are  L  sensors  that  measure  signals  impinging  on  the  array 
from  P  wideband  acoustic  sources  in  presence  of  spatially  uncorrelated  additive  noise.  The  sensor  array  can 
have  an  arbitrary  known  geometry.  In  this  part  of  the  report,  both  ULA  and  CSA  geometries  are  considered 
in  the  simulations.  The  sensor  measurement  is  denoted  by  Xi(n),  i  =  1, 2,  •  •  •  ,  L,  n  =  0, 1,  •  •  ■  ,  N  —  1, 
and  N  is  the  number  of  snapshots  that  will  be  used  in  the  DOA  estimation.  The  sources  are  denoted 
by  5i(n),  i  —  1,2,  ,P.  The  bandwidth  of  the  sources  is  divided  into  h  non-overlapping  subbands  or 

frequency  bins  of  interest  (the  way  to  select  the  h  subbands  or  frequency  bins  will  be  explained  on  later). 
Then,  at  each  frequency  bin  with  central  frequency  k  =  1, 2,  *  •  •  ,  h,  we  have  the  comp  lex- valued 
measurements 

XK)  =  [X1(a;fc))X2(wfe),.-.  ,XL(u,fc)],  k  —  1,2,---  ,h  (8) 
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The  narrowband  spatial  correlation  matrix  (SCM)  at  each  frequency  bin  is  then  computed  using 

RxM  =  X(uk)XH(uk),  *  =  1,2,...  (9) 


In  the  next  section,  we  present  a  brief  review  of  another  wideband  DOA  estimation  scheme  namely  IWM 
[33],  followed  by  the  WSF  algorithm  [51], [52].  The  steps  in  the  IWM  algorithm  may  be  summarized  as 
follows: 

Step  1:  Perform  eigenvalue  decomposition  (e.g.  Singular  Value  Decomposition  (SVD))  for  each  spatial 
correlation  matrix  to  obtain  the  noise  subspace  U (u^). 


Step  2:  For  a  bearing  angle  0,  evaluate 

P(uk,0)  = 


1 


where  a (a>fc,0)  is  the  steering  vector  that  depends  on  the  array  manifold. 

Step  3:  Average  the  power  over  all  the  selected  frequency  bins  for  the  bearing  angle  9 

,  h 


(10) 


(ii) 


k—l 


Step  4:  Calculate  P(@)  for  all  9  of  interest,  by  repeating  steps  2  and  3,  and  pick  those  for  which  P(-)  has 
a  peak  above  certain  level,  as  the  DOA  estimate  (s). 

4.3  Weighted  Subspace  Fitting  (WSF)  Algorithm 

The  steps  in  this  algorithm  can  be  summarized  as  follows: 

Step  1:  Perform  SVD  for  each  Rx(^fc)  to  obtain  the  respective  signal  subspace  Us(wfc).  Then  calculate 
the  diagonal  weighting  matrices  P(u>fc)  (rjk  x  rjk,  rjk: the  number  of  signal  components  contained  in 
R-x(^fc)  ) 


p(wfc)[fc,it] 


Afc(cUfc)  -  O-g 

V'MwfcW 


V  A/c  pik) 


(12) 


for  relatively  small  noise  variance  <j£.  X^s  are  the  eigenvalues  of  Rx(wjt).  If  the  subband  width  is 
zero  or  for  a  frequency  bin,  one  has  Tjk~  1  and  P(oJk)  w ill  be  a  scalar. 

Step  2:  Using  the  steering  vector  ©  =  [0i,  62,  •  *  *  ,  0p}T ,  compute 


C(0)  =  At(o;fe,©)Us(a;fc)P(a;fc) 

(13) 

where 

At(wfc,0)  =  {AH{uk,@)A(uk,&))  1AH(uk,@) 
A(uk,  ©)  =  [a(wfe,0i),a(wfe,02),--- ,a(wfc,0p)]T 

(14) 

(15) 

Step  3:  Calculate  the  following  criterion  for  all  the  steering  vectors  ©  of  interest,  by  repeating  the  steps 
1  and  2,  and  taking  the  one  that  minimizes  the  criterion  as  the  estimate  of  the  DOA  vector 

h 

@WSF  =  1™!^  |A(wfc,©)C(u>fc,©)  -  Us(wfc)P(u;fc) 
k—l 

2 

F 

(16) 
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4.4  Simulation  Results  of  Wideband  DOA  Estimation  Methods 

In  this  section  we  first  present  the  results  of  the  AR-based  method  for  generating  synthesized  data.  Figure 
4  shows  the  results  of  applying  IWM,  STCM  and  WSF  algorithms  using  a  uniform  linear  array  (ULA) 
to  a  synthetic  data  with  SNR  =  0  dB  generated  based  upon  the  truth  file  associated  with  neicrl502  in 
SAFE  II  data  set.  Figure  5  gives  the  results  of  the  three  DOA  estimation  methods  on  the  same  data 
but  using  a  circular  (wagon- wheel)  array  geometry  (CSA).  Several  intriguing  observations  can  be  made 
from  these  results.  Comparing  the  results  of  IWM  and  WSF  with  those  of  STCM  one  can  observe  that 
the  number  of  false  peaks  in  the  STCM  is  much  higher  than  the  other  two  algorithms.  Clearly,  WSF 
algorithm  is  the  superior  algorithm  as  it  provided  almost  perfect  DOA  estimation  and  tracking  ability  for 
both  array  configurations  (ULA  and  CSA)  and  in  this  very  low  SNR  condition.  Additionally,  the  results 
for  the  circular  array  has  much  smaller  mean  squared  error  (MSE)  as  compared  to  those  of  the  ULA.  This 
is  particularly  evident  by  comparing  Figures  4(c)  and  5(c)  for  WSF  algorithm. 

To  show  the  power  of  the  WSF  algorithm  for  resolving  closely  spaced  targets,  a  similar  experiment  was 
repeated  using  a  synthetically  generated  data  from  the  SAFE  II  data  set  associated  with  newrl507  with 
SNR  —  5  dB.  This  particular  case  presents  a  more  challenging  scenario  with  two  closely  spaced  targets 
that  move  in  abreast  formation.  Figures  6  and  7  show  the  results  of  IWM,  STCM  and  WSF  algorithms 
for  a  ULA  and  a  CSA,  respectively.  Once  again  the  results  indicate  the  superiority  of  the  WSF  method 
for  DOA  estimation  and  tracking  of  closely  spaced  targets.  To  demonstrate  the  real  usefulness  of  this 
algorithm  in  situations  where  more  than  two  abreast  targets  may  be  encountered,  we  repeated  the  same 
experiment  on  a  synthesized  data  set  with  three  closely  spaced  targets  in  an  abreast  formation  ( SNR  —  10 
dB).  Figures  8  and  9  show  the  results  of  these  algorithms  for  a  ULA  and  a  CSA  geometry,  respectively. 
These  results  clearly  demonstrate  the  superiority  of  WSF  algorithm  for  DOA  estimation  in  difficult  target 
formations  and  SNR  conditions.  We  can  also  conclude  that  the  superiority  of  the  WSF  algorithm  is  not 
dependent  on  the  number  of  the  targets  under  consideration,  as  we  have  consistently  achieved  very  good 
results  in  various  multiple  (>  2)  closely  spaced  target  scenarios. 


(a)  IWM 


(b)  STCM 


(c)  WSF 


Figure  4:  Uniform  Linear  Array  (ULA)  applied  to  synthetically  generated  SAFE  II  data  set  newrl502  with 
SNR  =0  dB.  The  x-axis  represents  time  in  s.  The  number  of  blocks  denotes  the  number  of  samples  plotted 
in  each  interval  of  5s. 


Next,  the  second  method  for  generating  realistic  synthesized  multiple  target  formations  is  applied. 
We  generated  three  cases  of  abreast,  single-file  and  staggered  targets.  In  all  these  cases,  the  radius  of 
the  array  was  4  ft.  and  the  center  was  located  at  (-5m, 5m).  For  the  abreast  case,  the  two  targets  are 
originally  at  (-300m, -300m)  and  (-300m, -350m).  For  the  single-file  case,  the  two  targets  are  originally  at 
(-300m, 100m)  and  (-250m, 100m).  In  the  staggered  formation  there  are  three  targets  at  initial  positions 
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Figure  5:  Cross  Shaped  Array  (CSA)  applied  to  synthetically  generated  SAFE  II  data  set  newrl502  with 
SNR  =0  dB.  The  z-axis  represents  time  in  s.  The  number  of  blocks  denotes  the  number  of  samples  plotted 
in  each  interval  of  5s. 


(a)  IWM  (b)  STCM  (c)  WSF 


Figure  6:  Uniform  Linear  Array  (ULA)  applied  to  synthetically  generated  SAFE  II  data  set  newrl507 
using  two  targets  with  SNR  =5  dB.  The  x-axis  represents  time  in  s.  The  number  of  blocks  denotes  the 
number  of  samples  plotted  in  each  interval  of  5s. 
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Figure  7:  Cross  Shaped  Array  (CSA)  applied  to  synthetically  generated  SAFE  II  data  set  newrl507  using 
two  targets  with  SNR  =5  dB.  The  z-axis  represents  time  in  s.  The  number  of  blocks  denotes  the  number 
of  samples  plotted  in  each  interval  of  5s. 


Figure  8:  Uniform  Linear  Array  (ULA)  applied  to  synthetically  generated  SAFE  II  data  set  newrl507 
using  three  targets  with  SNR  =10  dB.  The  z-axis  represents  time  in  s.  The  number  of  blocks  denotes  the 
number  of  samples  plotted  in  each  interval  of  5s. 
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Figure  9:  Cross  Shaped  Array  (CSA)  applied  to  synthetically  generated  SAFE  II  data  set  newrl507  using 
three  targets  with  SNR  =10  dB.  The  z-axis  represents  time  in  s.  The  number  of  blocks  denotes  the  number 
of  samples  plotted  in  each  interval  of  5s. 


(-300m, 100m),  (-250m, 100m),  (-275m, 150m).  In  all  these  cases,  the  targets  move  to  East  with  a  constant 
speed  of  25m/sec.  The  signal-to-noise  ratio  (SNR)  was  lOdB  for  these  cases,  which  represents  a  nominal 
operating  condition.  In  the  following,  only  the  results  of  the  STCM  and  WSF  algorithms  are  presented 
here  as  the  IWM  algorithm  did  not  produce  reliable  results. 

The  STCM  method  was  first  applied  to  these  cases.  The  results  are  shown  Figures  (10),  (11),  and  (12) 
respectively.  Figure  (10) (a)  shows  the  truth  tracks  of  the  two  abreast  targets.  Figure  (10) (b)  gives  the 
detected  DOA’s  for  this  case  using  the  STCM  method.  To  investigate  the  resolution  of  this  method  in 
separating  the  two  tracks,  we  generated  the  entire  response  of  the  STCM  before  peak  picking  operation. 
Figure  (10) (c)  shows  this  3-D  response  pattern  of  the  STCM  versus  time.  In  a  sense  this  plot  reflects  the 
beam  resolution  of  the  STCM  method.  The  frequency  range  used  for  the  detection  is  50 Hz  <  f  <  150 Hz. 
As  can  be  seen,  although  the  STCM  method  correctly  detects  the  track  of  the  two  targets,  it  could  not 
resolute  them  into  two  individual  targets.  Note  that  unlike  the  subband  MUSIC  algorithm  which  generates 
a  large  number  of  false  peaks,  the  STCM  method  does  not  have  this  problem  as  evident  from  both  Figures 
(10) (a)  and  (b).  Nonetheless,  if  we  expand  the  frequency  range  to  80 Hz  <  f  <  250 Hz,  we  start  seeing 
some  false  peaks  even  though  they  are  not  larger  than  those  of  the  targets  in  amplitude. 

Figures  (ll)(a)-(c)  show  the  results  of  the  single-file  case.  Figures  (12)(a)-(c),  on  the  other  hand,  show 
the  corresponding  plots  for  the  staggered  formation.  Again  as  we  can  see,  the  STCM  detects  the  track  of 
the  two  targets  but  could  not  resolute  them  well.  The  important  observation  of  this  study  is  that  the  STCM 
method  can  indeed  detect  the  DO  A  }s  of  groups  of  closely  spaced  targets .  This  implies  that  this  method  may 
be  a  very  useful  pre-detection  tool  that  provides  approximate  location  or  DOA  of  the  formation.  Once  this 
information  is  determined,  other  methods,  e.g.  WSF,  maybe  used  to  resolve  the  formation  into  several 
targets . 

In  contrast  to  the  STCM  that  is  unable  to  discriminate  between  the  number  of  targets,  the  WSF 
algorithm  is  capable  of  separating  the  targets  even  in  very  tight  formations.  To  demonstrate  this,  the  same 
synthetic  data  for  three  rather  difficult  scenarios,  namely  single-file,  staggered,  and  abreast  formations  are 
used.  The  results  shown  in  Figures  10,  11,  and  12  for  the  abreast,  single-file  and  staggered  formations, 
respectively  are  to  serve  as  reference  for  comparison  with  the  results  of  the  WSF  algorithm.  As  we  have 
stated  before,  Figures  10  and  11  showed  that  although  the  STCM  method  correctly  detects  the  track  of  the 
two  targets,  it  is  unable  to  resolute  them  into  two  individual  targets.  Similar  observation  can  be  drawn  from 
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Figure  10:  The  real  and  detected  tracks  of  two  abreast  targets. 


Figure  11:  The  real  and  detected  track  of  two  single-file  targets 
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Figure  12:  The  real  and  detected  track  of  three  staggered  targets 


the  results  of  the  staggered  case,  Figure  12,  with  three  closely  spaced  targets.  Since  the  STCM  algorithm 
is  computationally  more  efficient  comparing  to  the  WSF  algorithm  and  further  the  STCM  has  shown  good 
capabilities  in  identifying  the  target  tracks,  we  propose  to  use  the  STCM  algorithm  as  a  pre-processor  for 
the  WSF  algorithm.  In  this  way,  the  search  space  for  the  WSF  can  be  reduced  significantly,  yielding  a 
computationally  tractable  algorithm  with  very  good  resolution  capability  for  closely  spaced  targets.  This 
is  demonstrated  in  the  results  that  follow. 

Figures  13-15  depict  the  truth  tracks  of  the  abreast,  single-file  and  staggered  formations,  respectively 
together  with  the  detected  tracks  using  the  WSF  algorithm.  Note  that  the  actual  tracks  are  shown  by 
solid,  dotted  and  dashed  lines  in  these  figures  while  the  results  of  the  WSF  algorithm  are  shown  by  *,  x 
and  o  in  these  plots.  Comparing  these  results  with  those  of  the  STCM  method  clearly  show  the  advantages 
of  the  WSF  algorithm  in  resolving  closely  spaced  targets.  Among  future  research  topics  in  Phase  II  is  the 
design  and  development  of  a  computationally  efficient  algorithm  for  implementing  the  WSF  method. 


Figure  13:  The  WSF  algorithm  applied  to  the  synthetically  generated  Abreast  formation  data.  The  number 
of  blocks  in  the  x-axis  represents  time  in  s. 


Figure  14:  The  WSF  algorithm  applied  to  the  synthetically  generated  Single-File  formation  data.  The 
number  of  blocks  in  the  x-axis  represents  time  in  s. 
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Figure  15:  The  WSF  algorithm  applied  to  the  synthetically  generated  Staggered  formation  data.  The 
number  of  blocks  in  the  x-axis  represents  time  in  s. 


5  Array  Manifolds  &  DOA  Estimation 

In  this  part  of  the  study  we  concentrate  our  efforts  on  a  dual-size  array  geometry  that  consists  of  a 
rectangular  array  grid  to  be  spaced  much  greater  than  a  half-wavelength,  but  at  each  grid  point  a  5-element 
half- wavelength  spaced  cross- shaped  subarray  is  utilized  as  shown  in  Figure  16.  It  is  generally  known  that 
a  larger  array  aperture  will  be  able  to  produce  a  more  accurate  direction  of  arrival  (DOA)  estimate  as 
compared  to  a  smaller  array  aperture.  The  advantage  of  using  a  larger  array  aperture  is  also  manifested 
in  achieving  a  more  robust  resolution  of  very  closely  spaced  target  signatures.  In  order  to  formally  extend 
the  array  aperture  it  is  customary  to  employ  sparse  array  configurations  where  the  uniform  intra-array 
spacing  (denoted  by  A)  is  wider  than  the  half- wavelength  inter-subarray  spacing  (denoted  by  5  =  A/2). 
Therefore,  achieving  the  goal  of  more  accurate  estimate  without  using  extra  hardware.  Unfortunately, 
this  larger  intra-array  positioning  will  result  in  direction  cosines  u  and  v  estimates  (in  the  y  and  x  axes, 
respectively)  that  are  ambiguous.  To  resolve  this  ambiguity  a  sparse  but  regular  array  geometry  embedded 
with  a  spatial  invariance  may  be  employed  and  MUSIC-based  or  ESPRIT-based  DOA  algorithms  may  be 
developed  to  exploit  the  specific  geometry  and  to  resolve  the  ambiguities  [53]-[55]. 

The  smaller  spatial  invariance  6  =  A/2  would  result  in  an  unambiguous  but  high- variance  estimates  of 
the  DOAs  whereas  the  larger  spatial  invariance  A  >>  A/2  would  result  in  lower  invariance  at  the  expense 
of  cyclically  ambiguous  estimates  of  the  DOAs.  Consequently,  the  proposed  geometry  would  facilitate 
simultaneous  development  of  less  costly  hardware  architectures  and  computationally  efficient  algorithms 
that  would  take  advantage  of  each  spatial  invariance. 

For  the  sake  of  generality  and  completeness  we  consider  a  sparse  rectangular  array  geometry  as  an  L  x  M 
uniform  grid  consisting  of  five  (5)-element  identical  subarrays.  Consider  P  uncorrelated  and  narrowband 
signals  that  impinge  on  the  2-D  grid.  The  kth  signal  may  be  characterized  by  a  5 LM  x  1  array  manifold 


a(ufc,Ufc)  =  <n(uk,Vk)  <0  <is(uk,Vk),  k  =  l,---,P  (17) 

where  0  is  the  Kronecker  multiplication  operator,  q L,{ukivk)  is  a  LM  x  1  vector  consisting  of  the  spatial 
phase  factors  associated  with  the  sparse  uniform  array  grid,  qs(uk,vk)  is  a  5  x  1  vector  containing  the 
subarray  manifold,  and  uk  =  sin9kcos <pk>  vk  =  sinOkSinfik  represent  the  direction  cosines  along  the  y-axis 
and  x-axis,  respectively,  and  9k  and  4>k  represent  the  kth  source  elevation  (measured  wrt  the  vertical  2-axis) 
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Figure  16:  Extended  aperture  dual  invariance  sparse  array  with  rectangular  grid. 


and  azimuth  angles  (measured  wrt  the  positive  z-axis),  respectively.  For  the  rectangular  grid  geometry 
shown  in  Figure  16  we  have  the  following  relation 


e Quk 

ei2jLt-°Vk 

q  L(uk,Vk)  = 

[  VUk  \ 

ej^A(M-l)vk 

(18) 


where  Ax  »  A/2  and  Ay  »  A/2,  with  Ax  and  Ay  being  the  x-axis  and  y-axis  intra-array  distances 
between  the  grid  points  (nodes).  The  incident  signal  on  the  grid  may  now  be  represented  as 


z  (*)  = 

1 

_ 1 

=  ]a(ui,vi),---  ,a (up,vP)] 

'  «i(t)  ' 

_  zLM{t) 

:=A(ui,-  ,up,v i,*-*  >vp) 

sp(t ) 

(19) 


where  Zi(t),  i  =  1,  •  ■  •  ,LM  is  a  5  x  1  incident  signal  on  the  ith  subarray,  Si(t )  is  the  signal  associated 
with  source  i,  and  n(i)  refers  to  a  complex  valued  zero-mean  additive  white  noise.  Using  a  total  of  N 
snapshots  taken  at  distinct  times  {tn,n  -  l,---  ,N},  the  DOA  problem  is  then  reduced  to  determining 
{9k,<j>k, &  =  !,•••  ,  P}  from  the  5 LM  x  N  data  set  Z  =  [z(ti),  •  •  •  ,  z(tjy)]. 


5.1  Mathematical  Model  for  the  Subarray 

To  complete  the  MUSIC-based  algorithm  that  is  developed  here  for  the  half-wavelength  spaced  subarray 
at  each  grid  point  we  have 


q  s{uk,Vk) 


e~j(2n5x/\)vk  e—j(2nSy/X)uk  ,  ^  ejC2n5v /X)v,k ,  ej(2n5x /\)vk 


(20) 


and  where  6X  <  A/2  and  Sy  <  A/2  with  8X  and  Sy  denoting  the  x-axis  and  y-axis  distances  of  the  elements 
within  the  subarrays. 
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The  first  step  in  the  MUSIC-based  algorithm  is  to  compute  the  5 LM  x  P  signal-subspace  eigenvector 
matrix  Es  from  which  the  other  matrix  pencils  are  constructed  corresponding  to  the  P  largest  eigenvalues 
of  the  5 LM  x  5 LM  sample  spatial  correlation  matrix.  Namely,  we  have 


N 

Rzz  =  ZZH  =  J2  z(*")zH(*")  =  EsDsEf  +  EnDnE*  (21) 

71—1 

where  E„  is  the  5  LM  x  (5 LM  —  P )  noise-subspace  matrix  corresponding  to  the  5 LM  —  P  smallest  eigen¬ 
values,  Ds  is  the  P  x  P  diagonal  matrix  whose  diagonal  entries  are  the  largest  P  eigenvalues,  and  Dn  is  a 
(5 LM  —  P)  x  (5 LM  —  P )  diagonal  matrix  whose  diagonal  entries  are  the  5 LM  —  P  smallest  eigenvalues. 
As  the  number  of  snapshots  increase  the  value  for  Es  may  asymptotically  be  approximated  by 


Es  =  AT  =  [a(wi,i>i),  •  •  •  ,a.(up,vp)]T  (22) 

where  T  denotes  an  unknown  Px  P  nonsingular  matrix  to  be  determined.  Two  matrix  pencils  with  spatial 
invariance  are  now  constructed  along  the  y-axis  according  to 


A"  =  [/5(L— l)Mx5(L— 1)mI°5(L— 1)A1x5m]  A  (23) 

Ag  =  [05(L-1)Mx5mI^5(I-1)Mx5(L-1)m]  A  =  Ag(<l>“)t  (24) 

where  <f>,J  is  a  diagonal  matrix  with  P  diagonal  elements  and  f  is  the  complex  conjugate  operation.  The 
matrix  pencils  are  now  defined  as 

E1  =  [^5(l-l)Mx5(L-l)Ml05(L-l)Mx5M]  Es  (25) 

E2  =  [05(L-l)Mx5Ml^5(L-l)Mx5(L-l)M]  Es  (26) 

The  matrix  pencils  {Eg,  Eg}  are  formed  along  the  y-axis  by  having  all  (L  -  l)M  subarrays  at  the  top 
L  —  1  rows  on  the  overall  array’s  grid  to  contribute  data  to  Eg  and  by  having  all  (L  —  \)M  subarrays  at  the 
bottom  L- 1  rows  on  the  overall  array’s  grid  to  contribute  data  to  Eg.  The  above  two  matrices  are  related 
to  each  other  through  the  extended  spatial  invariance  Ay  »  A/2  and  can  yield  low  variance  but  cyclically 
ambiguous  DO  A  estimates  along  the  y-axis.  In  the  ideal  case  when  there  is  no  noise  or  the  number 
of  snapshots  is  sufficiently  large  we  may  write  Ej4 =  Eg,  so  that  <f*“  =  ((E'“)//E’1‘)_1((E,]‘);/Eg)  — 
(Tu)_1,5'uT“,  where  [^“]m  are  the  eigenvalues  of  <f>“. 

Prom  the  above  discussion  it  now  follows  that  a  set  of  highly  accurate  but  cyclically  ambiguous  estimates 
of  the  direction  of  cosines  along  the  y-axis  may  be  estimated  as  follows: 


Uj 


Z[*u]« 

2ttA2//A’ 


*  =  !,•••  ,P 


(27) 


where  Lz  refers  to  the  principal  argument  of  the  complex  number  z  in  the  range  of  —  tt  and  it.  Given  that 
Ay  »  A/2  and  -1  <ut  <1,  there  exists  a  set  of  cyclically  ambiguous  estimates  for  ui}  that  is 
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(28) 

(29) 


where  |Y|  and  [xj  denote  the  ceiling  and  floor  operators,  respectively.  Therefore,  T  may  be  estimated  to 
within  an  unknown  permutation  matrix  P“  as  T“  =  P“T. 

Similar  analysis  may  be  carried  out  for  the  extended  invariance  along  the  x-axis,  namely: 

=  ((Ei)hEi)-1((Ei)hE2)  =  (T")-1’®'1'^,  where  are  the  eigenvalues  of  and 
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The  pairs  {ui,Vi}  may  be  grouped  together  as  in  [53]-[55] . 


5.2  MUSIC-based  Ambiguity  Resolution 

In  order  to  overcome  and  remedy  the  cyclical  ambiguity  that  was  discussed  in  the  previous  section,  a 
MUSIC-based  approach  is  considered  in  this  part  of  the  project.  Specifically,  the  above  ambiguous  estimates 
of  (ui,Vi)  are  substituted  in  the  array  manifold  a(uk,vk).  The  best  candidate  is  now  selected  as  the  one 
that  minimizes  a  certain  distance  metric  between  E„  and  a(uk,vk).  This  disambiguation  takes  advantage 
of  the  superior  accuracy  of  the  MUSIC  algorithm,  but  avoiding  its  computational  complexity  and  cost. 
The  specifics  of  the  algorithm  are  as  follows: 


{uk,Vk}  =  arg  miri{UtV}\\E%(qL{u,v)®a(u,v))\\,  k  =  l,---,P  (33) 

The  simulation  results  reported  in  [53]-[55]  involving  two  closely  spaced  equal-power  narrowband  uncor¬ 
related  emitters  and  using  a  3  x  3  sparse  rectangular  arrays  have  shown  that  the  RMS  standard  deviation 
decreased  by  a  factor  of  about  50  folds.  Furthermore,  the  results  in  [54]  using  ESPRIT-based  approach  for 
a  9  x  9  sparse  grid  was  compared  with  a  9  x  9  half-wavelength  uniformly  spaced  square  array.  Although, 
in  terms  of  the  hardware  and  computational  complexities  the  two  array  geometries  are  comparable,  the 
DOA  estimation  performance  of  the  sparse  array  configuration  was  shown  to  be  surprisingly  better  than 
the  uniform  array,  which  also  failed  to  resolve  the  two  very  closely  spaced  sources. 

5.3  ESPRIT-based  DOA  Detection 

There  are  two  schemes  for  forming  the  subarrays  when  using  ESPRIT-based  method.  In  the  first  scheme 
we  let  all  the  sensors  at  the  same  position  within  the  cross  array  be  one  subarray,  e.g.  all  the  sensors  at 
the  center  of  the  cross  array  form  one  of  these  subarrays.  In  the  second  scheme  we  would  include  as  many 
sensors  into  one  subarray  as  possible  [53].  Specifically,  there  are  four  possible  subarrays  in  this  case:  (a) 
the  one  containing  the  top  two  sensors  on  each  vertical  leg  of  all  the  cross  arrays,  (b)  the  one  containing 
the  bottom  two  sensors  on  each  vertical  leg  of  all  the  cross  arrays,  (c)  the  one  containing  two  leftmost 
sensors  on  each  horizontal  leg  of  all  the  cross  arrays,  and  (d)  the  one  containing  two  right  most  sensors  on 
each  horizontal  leg  of  all  the  cross  arrays. 

The  advantages  of  the  first  ESPRIT  approach  over  the  MUSIC-based  approach  are  as  follows: 

1.  In  the  first  scheme,  there  are  only  ML  sensors  in  each  subarray,  thus  the  maximum  number  of  targets 
that  it  can  detect  is  ML  -  1,  while  in  the  second  scheme,  this  number  is  2 ML  -  1. 

2.  By  including  more  sensors  into  the  subarrays,  the  interference  from  noise  is  also  reduced. 
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3.  In  the  first  scheme,  there  are  5  subarrays.  Therefore,  integration  of  the  estimated  DOA’s  becomes 
an  issue,  while  in  the  second  scheme,  there  exists  an  effective  mechanism  described  in  [53], 

Due  to  the  above  reasons  in  our  implementation,  we  only  considered  the  ESPRIT-based  scheme. 

5.4  Simulation  Results 

Our  objective  here  is  to  further  investigate  the  use  of  the  above  techniques  for  development  of  ESPRIT- 
based  algorithms  that  have  specifically  taken  into  consideration  the  sparse  array  geometry  of  our  problem, 
and  to  benchmark  their  performance  against  the  other  standard  geometries. 

We  have  conducted  several  experiments  to  study  and  test  the  performance  of  the  above  method  using  the 
ESPRIT  approach.  The  setup  for  all  the  experiments  are  L  =  2,  M  —  2,  8X  —  5y  =  A/2,  Ax  =  Ay  =  10A. 
In  these  simulations,  we  assumed  the  targets  emit  narrowband  signals.  The  first  experiment  is  designed 
to  test  the  performance  of  the  detection  algorithm  as  a  function  of  the  variations  in  SNR.  The  result  is 
shown  in  Figure  17.  The  x-axis  measures  the  SNR  in  dBs  and  the  y-axis  denotes  the  Mean  Square  Error 
(MSE).  In  this  experiment,  two  targets  at  6\  —  62  =  tt/2,  <f>i  =  0.7n,  <j>2  =  0.97T  are  assumed.  From  the 

figure,  we  can  observe  that  as  expected,  the  MSE  decreases  when  the  SNR  increases. 


azumith  4>1=0.7tt:  <J>2  =  0.9rt 


Figure  17:  Performance  of  the  ESPRIT-based  algorithm  with  varying  SNR. 


The  second  experiment  is  designed  to  determine  the  DOA  detectability  of  this  method  for  resolving 
two  very  close  targets.  To  achieve  this,  the  DOA  of  the  first  target  was  fixed,  and  the  DOA  of  the  second 
target  was  continuously  increased  by  0.017T  increments.  The  results  in  Figure  18  reveals  an  interesting 
behavior.  That  is,  with  an  exception  of  several  specific  DOA  differences  (specifically  at  0.047r,0.097r,  etc.), 
this  method  exhibits  very  good  performance  in  detecting  two  closely  spaced  targets. 

To  test  the  capability  of  the  algorithm  for  detecting  multiple  (>  2)  closely  spaced  targets,  in  the 
next  experiment,  a  4  x  4  grid  was  used.  There  are  four  targets  with  the  corresponding  DOA’s  that  are 
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0.l7r,0.27r,0.37r,0.47r,  respectively.  The  results  were  obtained  by  rounding  the  detected  angle  to  multiples 
of  0.1-7T.  From  the  results  presented  in  Figure  19,  it  can  easily  be  observed  that  this  method  can  detect  the 
true  DOA’s  very  frequently.  Future  work  should  include  extending  this  approach  to  wideband  sources. 

6  Source  Separation  Methods 

In  this  section  we  investigate  various  ICA  algorithms  and  their  applicability  for  blind  source  separation  of 
multiple  targets  moving  in  closely  spaced  staggered  and  abreast  formations.  Among  the  researched  algo¬ 
rithms  are:  instantaneous  ICA  (IICA)  [56]-[61]  and  convolutive  (CICA)  [62]-[65].  We  will  also  investigate 
neural  network-based  ICA  methods,  namely  fast  ICA  in  [66],  The  actual  data  from  SAFE  II  database 
as  well  as  the  synthesized  data  containing  from  2-5  closely  spaced  targets  are  utilized  to  determine  the 
performance  of  the  IICA-based  algorithm.  It  was  observed  that  when  the  number  of  the  targets  for  a  given 
node  is  lower  than  the  number  of  the  microphones,  the  IICA  algorithm  may  not  be  able  to  always  correctly 
“de-confuse”  or  “discriminate”  among  the  targets  in  certain  situations.  Consequently,  our  future  research 
would  involve  modification  or  enhancement  to  the  existing  IICA  algorithm  to  remedy  and  eliminate  this 
“ambiguity”  in  determining  the  number  of  the  targets.  Furthermore,  we  shall  investigate  the  applicability 
of  CICA-based  algorithms  to  the  real  (SAFE  II  data)  as  well  as  the  synthetic  data  in  order  to  determine 
the  advantage  and  disadvantages  of  the  two  methods  in  terms  of  their  “de-confusing”  and  “discriminatory” 
ability  in  resolving  multiple  closely  spaced  targets. 

We  first  introduce  one  of  the  well-known  instantaneous  ICA  (IICA)  algorithms  based  on  the  natural 
gradient  update  of  the  demixing  matrix.  A  number  of  nonlinear  estimating  functions  that  affect  the 
performance  of  the  algorithm  are  also  given.  Simulation  results  are  then  provided  to  demonstrate  the 
performance  of  this  algorithm.  Finally,  we  present  a  convolutive  ICA  (CICA)  algorithm  that  is  obtained  by 
extending  the  above-mentioned  IICA,  and  also  mention  some  other  promising  CICA  algorithms  developed 
in  the  literature. 

6.1  Natural  Gradient  Based  IICA  Algorithm 

The  goal  of  ICA  or  blind  source  separation  (BSS)  is  to  estimate  P  statistically-independent  source  signals 
from  L  different  measured  mixtures  of  these  signals.  The  linearly  mixed  measurements  are  obtained  by  a 
passive  acoustic  sensor  array  with  L  sensors,  where  P  <  L.  For  instantaneous  mixing,  a  measured  signal 
vector  x(k)  —  [2q(fc),  •  •  •  ,  xi(k)]T  is  assumed  to  be  linearly  mixed  from  the  unknown  source  signal  vector 
s(k)  =  [si(fc),  •  •  •  ,  sp(k)]T,  that  is 

x(k)  =  As  (k)  (34) 

where  A  is  an  unknown  ( L  x  P)  mixing  matrix.  The  most  important  assumption  is  the  source  indepen¬ 
dence;  i.e.,  all  Si(k )  are  statistically-independent  of  each  other.  For  the  linear  mixing,  intuitively  a  linear 
demixing  model  may  be  applied  for  the  source  separation  task,  that  is 

y (*)  =  W (fc)x(fc),  fc  =  0,l,---  ,N-1  (35) 

y(*0  =  [yi{k),y2(k),---  ,yp{k)]T 

where  W(fc)  is  a  P  x  L  demixing  matrix  with  adjustable  elements  Wij(k).  We  now  seek  for  an  algorithm 
to  determine  the  demixing  matrix  W (k)  such  that  the  combined  matrix  C(fc)  =  W(fc)A  is  mapped  to 
3>D,  where  is  an  P  x  P  permutation  matrix  and  D  is  a  diagonal  nonsingular  scaling  matrix.  Clearly, 
the  output  of  the  demixer  yi(k)  will  be  one  of  the  scaled  source  signals.  Ordering  and/or  identifying  the 
recovered  source  signals  y(k)  is  not  possible  without  additional  information  about  them. 

Numerous  techniques  have  been  proposed  to  solve  this  problem  (see  [56]- [61]).  A  density-based  ICA 
method  known  as  natural  gradient  algorithm  proposed  by  Amari  [58]  is  of  particular  interest  owing  to  its 
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good  performance,  robust  stability  properties,  and  low  computational  requirements  in  comparison  with 


other  similar  algorithms.  Specifically,  it  is  given  by 

W(fc  +  1)  =  W(fc)  +  fi{k)  [I  -  f  (y (k))yT (k)]  W(fc)  (36) 

where  g(k)  is  a  positive  learning  rate  parameter,  f(-)  is  a  vector  of  a  scalar  nonlinear  function  /(•)  known 
as  estimating  function  or  scoring  function.  Several  possible  candidates  for  this  function  are  listed  below. 

f{y)  —  sgn(y )  ( instantaneously  —  mixed  speech  sources )  (37) 

f(y)  —  (3y  +  y3,  (3  >  0  ( negative  —  kurtosis  sources )  (38) 

f(y)  =  (3y  +  tanh(m/),  /3  >  0  ( super  —  Gaussian  sources)  (39) 

f(y )  =  y y 3  ~  Ty5  “  Yy7  +  Ty9  +  \yl1  ( 9eneral )  (40) 

f{y)  =  y5  +  | y7  +  y?/9  +  y2/U  -  yV3  +  128^15  “  V17  {general)  (41) 


We  have  found  that  when  applied  to  the  above  algorithm  (36),  the  estimating  function  of  higher  order 
requires  more  snapshots  or  data  to  ensure  convergence  and  stability.  Generally,  it  is  preferable  to  use  a 
low-order  function  as  long  as  the  selected  function  has  a  good  fit  to  the  statistical  properties  of  the  sources. 

Figure  20  shows  the  three  independent  sources  and  their  estimates  by  the  IICA  algorithm  (36)  with 
estimating  function  (41).  Clearly,  the  sources  are  well-separated  after  300  batch  iterations.  In  Figure  21, 
the  rejection  indices  of  the  algorithm  using  three  different  estimating  functions  are  also  shown.  It  can  be 
observed  that  the  estimating  functions  (40)  and  (41)  lead  to  the  same  convergence  rate,  but  the  latter 
yields  smaller  steady-state  rejection  index  on  average.  However,  (41)  is  more  complicated  and  needs  more 
computational  time.  Estimating  function  (38)  is  much  simpler  than  (40)  and  (41).  However,  it  takes  longer 
(by  about  50  iterations)  to  converge.  Note  that  the  use  of  the  estimating  functions  (37)  and  (39)  was  not 
successful  in  our  simulations  for  this  example.  Figure  22  shows  one  example  of  using  the  function  (39). 
The  corresponding  convergence  of  the  rejection  index  is  given  in  Figure  21  (d),  where  the  learning  rate 
parameter  was  adjusted  to  0.075  such  that  the  convergence  is  similar  to  those  using  estimating  functions 
(40)  and  (41).  Other  parameters  /?  and  a  were  set  to  0.1  and  0.5,  respectively.  Adjusting  these  two 
parameters  did  not  change  the  performance  of  the  algorithm  significantly.  Estimating  function  (37)  just 
simply  failed  in  all  the  runs. 

6.2  Convolutive  ICA  (CICA)  Algorithm 

In  many  applications,  the  use  of  convolutive  mixing  is  more  practical,  since  sources  impinge  on  each  sensor 
at  different  times.  Clearly,  if  the  source  frequencies  and  the  sensor  spacings  are  much  smaller  than  the 
travelling  speed  of  source,  the  time  lag  or  delay  may  be  ignored,  and  the  IICA  may  be  approximately 
applied  equally  effectively.  On  the  other  hand,  if  this  is  not  the  case,  we  may  need  to  consider  using  CICA- 
type  algorithms.  The  convolutive  mixing  is  more  complicated  than  the  instantaneous  one.  In  convolutive 
mixing,  we  have 


oo 

x(k)  =  Ajs(k  -  j)  (42) 

j—  oo 

where  the  sequence  of  matrices  A j  ( L  x  P)  is  the  multichannel  impulse  response  of  the  unknown  mixing 
system.  The  demixing  system  is  assumed  to  be  governed  by 

M 

y(k)  =  £W j(k)x(k-j)  (43) 

j=o 
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Figure  20:  Three  actual  sources  used  (left  column)  and  their  estimates  (right  column)  by  algorithm  (36) 
and  estimating  function  (41). 


The  CICA  algorithm  is  obtained  by  extending  the  spatial-only  scheme  of  (36)  to  the  spatio-temporal 
case  and  by  introducing  certain  approximations,  which  may  then  be  represented  by 

W j(k  + 1)  =  W j(k)  +  fj.(k)  [Wj(fc)  -  f(y(fe  -  M))uT(k  -  j)]  (44) 

M 

n{k)  =  YJ^li-q{k)y{k-q)  (45) 

9= 0 

The  above  algorithm  has  an  elegant  form  and  is  computationally  efficient.  However,  some  of  the 
approximations  made  in  order  to  derive  it  may  degrade  the  performance.  This  is  an  area  that  needs  much 
further  research.  Recently,  several  other  promising  techniques  have  been  proposed  in  the  literature.  The 
geometric  source  separation  method  [62],  [63]  and  the  generalized  sidelobe  decorelator  [62],  [64]  based 
on  the  cross-power  spectral  minimization  can  resolve  the  ambiguities  inherent  in  the  CICA  algorithm 
by  introducing  geometric  constraints.  A  new  CICA  method  has  also  been  proposed  in  [65],  which  uses  a 
forward  FIR  model  to  identify  the  multi-path  channel  and  a  backward  FIR  one  to  generate  the  independent 
sources.  All  these  methods  have  been  confirmed  effective  when  applied  to  acoustic  signals  recorded  in  a 
reverberant  environment.  Their  applicability  and  suitability  to  our  problem  is  under  investigation. 

6.3  Simulation  Results  of  IICA  Algorithm  on  Actual  SAFE  II  Data 

Our  preliminary  results  of  the  IICA  algorithm  on  two  sets  of  real  SAFE  II  data  are  given  below.  In  all  the 
subsequent  figures  the  plots  on  the  left  column  correspond  to  the  sensor  array  measurements  x*  (except 
the  last  row),  and  the  ones  on  the  right  column  designate  the  demixer  output  yi  (except  the  last  row).  The 
plots  in  the  last  row  show  the  convergence  of  the  weights.  Another  performance  indicator  is  the  mutual 
(1st  order)  correlation  coefficient  matrix  (CCM)  evaluated  at  steady-state  condition  of  the  algorithm.  The 
learning  rate  is  set  to  be  0.075  and  the  number  of  snapshots  is  2000. 

The  results  shown  in  Figures  23  to  28  are  applied  to  the  actual  data  for  Node  1  (newr!502  in  SAFE 
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Figure  22:  Three  actual  sources  used  (left  column)  and  their  estimates  (right  column)  by  algorithm  (36) 
and  estimating  function  (39). 


II)  containing  two  targets  using  the  estimating  functions  (37)-(40).  The  corresponding  CCM’s  for  these 
functions  are  as  follows,  respectively 
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-0.040 
1.000  ’ 
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1.000 
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The  results  for  the  estimating  function  (41)  are  not  shown  as  the  algorithm  did  not  converge  for  this 
choice  on  the  real  SAFE  II  data.  This  could  be  due  to  low  number  of  snapshots  used  during  the  training 
process.  From  the  above  results  it  can  be  observed  that  the  algorithm  using  functions  (37)  and  (39)  is 
capable  of  identifying  and  separating  the  two  targets.  However,  the  performance  is  not  as  satisfactory 
using  functions  (38)  and  (40).  These  observations  follow  by  noting  the  very  small  off-diagonal  elements 
for  the  2  sensors  cases  for  (37)  and  (39),  and  relatively  significant  off-diagonal  elements  for  the  cases  (38) 
and  (40).  On  the  other  hand,  separating  the  array  signals  into  three  distinct  signals  results  in  a  relatively 
large  off-diagonal  elements  in  the  3  sensor  case  as  seen  from  COM  matrices  for  the  functions  (37),  (38), 
and  (40) .  This  may  be  viewed  as  an  indication  that  the  actual  number  of  sources  is  two  in  these  cases  as 
opposed  to  three.  However,  this  conclusion  is  weakened  by  noting  that  the  COM  in  case  of  3  sensors  for 
(39)  suggests  the  presence  of  three  targets  contrary  to  the  reality  in  this  case.  This  clearly  indicates  that 
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further  research  is  warranted  to  resolve  this  ambiguity  in  determining  the  true  number  of  targets  in  the 
collected  sensory  data. 

The  results  tested  for  the  actual  data  for  Node  1  (newrl536)  from  SAFE  II  containing  five  targets  and 
using  two  to  five  sensors  are  shown  in  Figures  29  to  32,  respectively.  The  CCM’s  for  various  choices  of 
sensors  and  using  the  estimating  function  (38)  are  shown  below.  It  can  be  observed  that  the  algorithm  is 
capable  of  identifying  and  perhaps  separating  the  targets  as  indicated  by  the  very  small  off-diagonal  terms 
in  the  CCM’s.  The  results  for  the  other  functions  are  not  as  good  as  these  and  in  some  cases  the  weights 
did  not  converge  at  all. 
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Several  interesting  observations  can  be  made  from  these  preliminary  results.  These  are  given  below. 

1.  When  the  number  of  independent  sources,  P,  is  larger  than  the  number  of  sensors,  L,  i.e.  P  >  L, 
the  IICA  algorithm  under-estimates  the  number  of  sources  and  identifies  a  maximum  of  L  sources. 
This  is  due  to  the  fact  that  the  algorithm  may  treat  a  particular  grouping  or  lumping  of  independent 
sources  as  a  single  independent  source.  On  the  other  hand,  when  L>  P,  the  algorithm  over-estimates 
the  number  of  sources.  In  this  case,  in  addition  to  the  independent  sources,  linear  combination  of 
them  may  also  be  identified  as  potential  sources.  It  is  indeed  possible  to  narrow  down  the  choices 
by  performing  some  type  of  postprocessing  on  the  extracted  sources.  This  is  one  aspect  that  would 
require  further  research  and  modifications  for  a  successful  implementation  of  the  algorithm. 

2.  An  objective  performance  evaluation  of  the  results  of  this  algorithm  seems  to  be  rather  difficult,  since 
the  outcome  of  IICA  is  highly  affected  by  several  factors  including:  nature  of  the  sources,  choice  of 
the  estimating  function  used,  and  actual  mixing  situation  at  the  array,  i.e.,  instantaneous  mixing  or 
convolutive  mixing.  One  of  our  objectives  and  goals  for  future  (Phase  II)  is  to  conduct  further  and 
more  thorough  investigation  of  this  scheme  and  provide  a  comprehensive  benchmark  with  the  other 
source  separation  and  DOA  estimation  algorithms.  Of  great  importance  to  this  study  is  to  determine 
how  this  method  performs  when  scenarios  with  tight  target  formations  are  encountered. 

6.4  Fast  ICA  Algorithm 

In  this  part  of  the  report  we  investigate,  develop  and  test  a  Fast  ICA-based  algorithm  [66] .  Fast  ICA 
algorithm  is  an  algorithm  that  attempts  to  maximize  the  so-called  contrast  function  (neg-entropy)  required 
to  implement  ICA  [66],  Below  we  follow  a  step-by-step  process  for  developing  and  implementing  this 
algorithm.  This  algorithm  is  then  applied  the  synthetically  generated  data  set  as  well  as  to  the  real  SAFE 
II  sets  corresponding  to  series  503,  510  and  537  scenarios.  Interesting  results  are  obtained.  Also  other 
interesting  research  venues  are  identified  and  indicated.  One  aspect  that  we  intend  to  investigate  in  Phase 
II  will  be  the  development  of  data  association  and  fusion  schemes  that  attempt  to  reduce  the  ambiguities 
that  may  arise  from  the  application  of  Fast  ICA  algorithm  to  different  nodes  in  the  array  configuration. 
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Figure  23:  Two  sensor  array  signals  and  two  demixer  outputs  case  using  estimating  function  (37)  for  SAFE 
II  real  data  newr!502. 
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Figure  24:  Three  sensor  array  signals  and  three  demixer  outputs  case  using  estimating  function  (37)  for 
SAFE  II  real  data  newrl502. 
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Figure  25:  Two  sensor  array  signals  and  two  demixer  outputs  case  using  estimating  function  (38)  for  SAFE 
II  real  data  newr!502. 
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Figure  26:  Three  sensor  array  signals  and  three  demixer  outputs  case  using  estimating  function  (38)  for 
SAFE  II  real  data  newr!502. 
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Figure  27:  Two  sensor  array  signals  and  two  demixer  outputs  case  using  estimating  function  (39)  for  SAFE 
II  real  data  newr!502. 
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Figure  28:  Three  sensor  array  signals  and  three  demixer  outputs  case  using  estimating  function  (39)  for 
SAFE  II  real  data  newr!502. 
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Figure  29:  Two  sensor  array  signals  and  two  demixer  outputs  case  using  estimating  function  (38)  for  SAFE 
II  real  data  newr!536. 
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Figure  30:  Three  sensor  array  signals  and  three  demixer  outputs  case  using  estimating  function  (38)  for 
SAFE  II  real  data  newr!536. 
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Figure  31:  Four  sensor  array  signals  and  Four  demixer  outputs  case  using  estimating  function  (38)  for 
SAFE  II  real  data  newrl536. 
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Figure  32:  Five  sensor  array  signals  and  Five  demixer  outputs  case  using  estimating  function  (38)  for 
SAFE  II  real  data  newrl536. 
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As  mentioned  before,  the  estimation  of  data  model  of  ICA  is  usually  performed  by  formulating  an 
objective  function  and  then  minimizing  or  maximizing  it.  Often  such  a  function  is  called  a  contrast 
function.  An  ICA  method  can  be  expressed  as  the  combination  of  an  objective  function  and  its  optimization 
algorithm.  The  statistical  properties  of  the  ICA  method  depend  on  the  choice  of  the  objective  function, 
and  the  algorithmic  properties  depend  on  the  optimization  algorithm.  These  two  classes  of  properties  are 
independent  in  the  sense  that  different  optimization  methods  can  be  used  to  optimize  a  single  objective 
function,  and  conversely  a  single  optimization  method  may  be  used  to  optimize  different  objective  functions 
[61].  The  key  to  estimating  the  ICA  model  is  non-Gaussian  property.  Two  classes  of  contrast  functions 
are  available  to  estimate  the  data  model  in  ICA,  which  are  described  below. 

A:  Multi-Unit  Contrast  Functions 

Using  these  type  of  contrast  functions  the  whole  data  model  of  ICA  is  estimated  in  one  shot. 

•  Likelihood  and  network  entropy:  The  likelihood  in  the  noise-free  ICA  model  is  formulated  and  then  the 
model  is  estimated  by  a  maximum  likelihood  (ML)  method  [67],  From  a  neural  network  viewpoint 
this  is  equivalent  to  maximizing  the  output  entropy  (information  flow)  of  a  neural  network  with 
nonlinear  outputs.  Often  this  is  termed  as  “infomax”  [57], [60]. 

•  Mutual  information  and  Kullback-Leibler  divergence:  The  mutual  information  is  a  natural  measure 
of  the  dependence  between  random  variables.  It  is  always  non-negative,  and  zero  if  and  only  if  the 
variables  are  statistically  independent.  Finding  a  transform  that  minimizes  the  mutual  information 
between  the  components  Si  is  a  very  natural  way  of  estimating  the  ICA  model  [56]. 

Among  other  contrast  functions  considered  in  the  literature  are  nonlinear  cross  correlation  by  Jutten 
and  Herault  [68],  nonlinear  PC  A  criteria  [69],  higher-order  cumulant  tensors  [70]  and  weighted  covariance 
matrix  [71], 

B:  One-Unit  Contrast  Functions 

These  functions  estimate  single  component  at  a  time  instead  of  the  whole  data  model.  This  approach 
shows  clearly  the  connection  to  neural  networks  and  prior  knowledge  of  the  number  of  independent  com¬ 
ponents  is  not  needed.  After  estimating  one  independent  component,  one  can  use  simple  decorrelation  to 
find  a  different  independent  component,  since  the  independent  components  are  by  definition  uncorrelated. 

Negentropy:  A  most  natural  information-theoretic  one-unit  contrast  function  is  negentropy.  The  indepen¬ 
dent  components  correspond  to  directions  in  which  the  differential  entropy  of  wTx  is  minimized.  However, 
a  modification  has  to  be  made,  since  differential  entropy  is  not  invariant  for  scale  transformations.  In  order 
to  obtain  a  linearly  invariant  version  of  entropy,  the  negentropy  J  is  defined  as 

J(y)  =  H(ygauss )  -  H{ y)  (46) 

where  y gauss  is  a  Gaussian  random  vector  of  the  same  covariance  matrix  as  y  and  H  denotes  the  entropy. 
Negentropy  is  always  negative,  and  zero  if  and  only  if  y  has  a  Gaussian  distribution  [56].  The  usefulness 
of  this  definition  is  prominent  when  mutual  information  is  expressed  using  negentropy  as 

1  TlCy 

/(».,'  •  • ,  yp)  =  J(y)  -  E  •’<■*>  +  2 log  <47) 

i 

where  Cy  is  the  covariance  matrix  of  y,  and  the  Cyt  are  its  diagonal  elements.  If  the  iji's  are  uncorrelated, 
the  third  term  is  0,  and  hence  we  get 

i  =  J(y)  ~  J(Vi)  (48) 
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As  negentropy  is  invariant  for  linear  transformation,  it  is  obvious  that  finding  maximum  negentropy  di¬ 
rections,  i.e.,  directions  where  the  elements  of  the  sum  J(y;)  are  maximized,  is  equivalent  to  finding  a 
representation  in  which  mutual  information  is  minimized  [56],  The  estimation  of  negentropy  is  difficult, 
and  therefore  this  contrast  function  remains  a  topic  for  theoretical  investigation.  Negentropy  can  be  ap¬ 
proximated  by  higher-order  cumulants  or  by  more  accurate  methods  as  shown  below. 

As  in  the  multi- unit  case,  negentropy  can  be  approximated  by  higher-order  cumulants  such  as  [72]: 

J( y)  «  ^«s(y)2  +  -^(y)2  (49) 

where  kx  is  the  z-th  order  cumulant  of  y.  However,  the  validity  of  such  approximation  may  be  rather 
limited.  In  [72],  it  was  argued  that  cumulant- based  approximations  of  negentropy  are  inaccurate,  and  in 
many  cases  too  sensitive  to  outliers.  The  new  approximations  are  based  on  the  maximum- entropy  principle. 
In  general,  one  may  obtain  an  approximation  of  the  form  J(y)  ~  X)f=i  ci[E{Gi(y)}  —  E{Gl(i>)}]2  where 
Cj’s  are  some  positive  constants,  v  is  a  Gaussian  variable  of  zero  mean  and  unit  variance,  and  G,’s  are 
practically  any  non-quadratic  functions.  In  the  case  where  one  uses  only  one  function  G,  the  approximation 
becomes  of  the  form: 

J(y)*c[E{G(y)}-E{G(v)}}2  (50) 

The  result  is  a  generalization  of  the  cumulant-based  approximation  of  (51)  where  y  is  symmetric.  For 
instance,  using  G(y)  =  y4,  one  then  obtains  (51)  from  (52).  This  one-unit  contrast  function  will  be  used 
later  for  the  development  of  Fast  ICA  algorithm. 

6.4.1  Preprocessing  for  ICA 

Before  applying  an  ICA  algorithm  on  the  observed  data,  it  is  usually  very  useful  to  do  some  preprocessing. 
The  steps  involved  are  given  below. 

1.  Centering-.  The  mean  vector  of  x,  denoted  by  m  =  £[x],  is  subtracted  from  x  to  make  it  a  zero-mean 
variable.  This  preprocessing  is  done,  to  simplify  the  ICA  algorithm.  After  estimating  the  mixing 
matrix  A  with  centered  data,  the  estimation  can  be  completed  by  adding  the  mean  vector  of  s  back 
to  the  centered  estimates  of  s.  The  mean  vector  of  s  is  given  by  A-1m. 

2.  Whitening:  After  centering  the  observed  vector  x,  it  is  transformed  linearly  so  that  the  new  vector 
is  white  i.e.,  its  components  are  uncorrelated  and  their  variances  equal  to  unity.  In  other  words,  the 
covariance  matrix  of  x  equals  the  identity  matrix: 

£[yyr]  =  I  (51) 

The  popular  method  of  whitening  is  to  use  the  eigenvalue  decomposition  (EVD)  of  the  covariance 

matrix  E[x.xT}  =  HDHT,  where  H  is  the  orthogonal  matrix  of  eigenvectors  of  F[xxT]  and  D  is  the 
diagonal  matrix  of  its  eigenvalues,  D  =  Diag{d\ ,  •  •  •  ,  dp).  Now 

x  =  EH~^Etx  (52) 

Whitening  reduces  the  number  of  parameters  to  be  estimated  and  it  transforms  the  mixing  matrix 
into  a  new  one,  A.  From  equations  (34)  and  (52) 

x  =  EKHetAs  =  As  (53) 
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6.4.2  Fast  ICA  for  One  Unit 


By  “unit”  we  imply  an  artificial  neuron  having  a  weight  vector  w  where  the  neuron  is  allowed  to  update  its 
weight  through  a  learning  rule.  The  Fast  ICA  learning  rule  finds  a  direction,  i.e.,  a  unit  vector  w  such  that 
the  projection  wrx  maximizes  nongaussianity.  Nongaussianity  here  is  measured  by  the  approximation  of 
negentropy  J(wTx)  given  as  J(y)  =  c{E[G( y)]  -  E[G(is)}}2.  G  is  practically  any  non-quadratic  function, 
c  is  any  The  algorithm  is  now  based  on  a  fixed-point  iterative  scheme  for  finding  a  maximum  of  the 
nongaussianity  of  wTx.  Let  g  denote  the  derivative  of  the  nonquadratic  function  G  used  above.  The  basic 
form  of  the  Fast  ICA  algorithm  can  be  expressed  as  follows  [66]: 

1.  Choose  an  initial  (e.g.  random)  weight  vector  w(0), 

2.  Define  w+  =  E[x.g( wTx)]  -  F[x#(wrx)]w,  and 

3.  Let  w  =  w+/||w+|[,  denote  the  normalized  w+ 

4.  If  there  is  no  convergence,  then  go  back  to  step  2. 

Here  by  convergence  we  imply  that  the  old  and  the  new  values  of  w  point  in  the  same  direction,  i.e. 
their  dot-product  is  (or  almost)  equal  to  1.  the  same  In  the  derivation  of  the  Fast  ICA,  first  the  maxima 
of  the  approximation  of  the  negentropy  of  wrx  are  obtained  at  certain  optima  of  E[G( wTx)].  According 
to  the  Kuhn-Tucker  conditions,  the  optima  of  E[G(wTx)]  under  the  constraint  £[(wTx)2]  =  ||w||2  =  1  are 
obtained  at  points  where 

l?[x<7(wTx)}  —  /?w  —  0  (54) 


Denoting  the  function  on  the  left  side  of  (54)  by  F,  its  Jacobian  matrix  becomes 

JF( w)  =  F[xxT£7(wTx)]  -  j3I 

Since  the  data  is  whitened,  a  reasonable  approximation  seems  to  be 

E[xxTgf(  wTx)]  =  F[xxT]F[#'(wrx)]  =  E[g'(wTx.)]I 


(55) 


(56) 


Thus,  the  Jacobian  matrix  becomes  diagonal  and  can  easily  be  inverted.  The  following  Newton’s 
iteration  is  approximated: 

w+  =  w  -  |F[x^(wrx)]  -  {3w\/\E[g'(wTx)\  -  (3\  (57) 

This  again  can  be  simplified  by  multiplying  both  sides  of  (57)  by  /3  -  E[g/(wTx)\.  This  gives  after  some 
algebraic  manipulation,  the  Fast  ICA  iteration.  In  practice,  the  expectations  in  Fast  ICA  must  be  replaced 
by  their  estimates.  The  natural  estimates  are  clearly  the  corresponding  sample  means.  The  averages  can 
be  estimated  using  a  smaller  sample,  whose  size  may  have  a  considerable  effect  on  the  accuracy  of  the  final 
estimates.  The  sample  points  should  be  chosen  separately  at  every  iteration.  If  the  convergence  is  not 
satisfactory,  an  increment  of  the  sample  size  may  be  required. 
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6.4.3  Fast  ICA  for  Several  Units 

To  estimate  several  independent  components  the  one-unit  Fast  ICA  need  to  be  run  using  several  units  with 
weight  vectors  wi,  •  •  •  ,  w p.  To  prevent  different  vectors  from  converging  to  the  same  maxima  the  outputs 
wf  x,  •  •  •  ,  WpX  must  be  decorrelated  after  every  iteration.  A  simple  way  of  achieving  decorrelation  is  a 
deflation  scheme  based  on  a  Gram-Schmidt  like  decorrelation.  After  estimating  q  independent  components, 
or  q  vectors  wi,  •  •  •  ,wg,  the  one-unit  fixed-point  algorithm  is  run  for  wp+i  and  after  every  iteration  step 
the  “projections”  wLwjWj,  j  —  1,  •  •  •  ,  q  of  the  previously  estimated  q  vectors  are  subtracted  from  w9+i, 
and  the  result  is  then  normalized  as  follows: 

1.  Define  wg+i  =  wg+i  -  wJ+iwjCwj,  and  then  normalize  it  according  to 

2.  w9+i  =  Wg+i/fw^Cw^i]1/2 

where  C  is  the  Covariance  of  the  data.  Since  the  data  is  sphered,  therefore  C  is  replaced  by  an  identity 
matrix  I. 

Properties  of  Fast  ICA  Algorithm 

The  Fast  ICA  algorithm  and  its  contrast  functions  have  a  number  of  desirable  properties: 

•  The  convergence  is  cubic  or  at  least  quadratic.  This  is  in  contrast  to  other  ICA-based  schemes  that 
are  derived  based  on  stochastic  gradient  descent  methods,  where  the  convergence  is  only  linear.  This 
means  a  very  fast  convergence. 

•  The  algorithm  is  easy  to  use  as  there  is  no  step  size  to  choose. 

•  It  finds  directly  independent  components  of  any  non-Gaussian  distribution  using  any  nonlinearity. 

•  The  performance  of  the  method  can  be  optimized  by  choosing  a  suitable  nonlinearity. 

•  The  independent  components  can  be  estimated  one-by-one.  This  is  useful  in  exploratory  data  analy¬ 
sis,  and  decreases  the  computational  load  of  the  method  in  cases  where  only  some  of  the  independent 
components  need  to  be  estimated. 

•  The  Fast  ICA  algorithm  enjoys  many  advantages  of  neural  networks  algorithms  such  as:  parallel  and 
distributed  processing,  simplicity,  and  requiring  little  memory  space. 

6.5  Simulation  Results  of  Fast  ICA  on  Synthesized  Data 

The  ICA  is  performed  for  each  window  of  1  second  duration.  We  have  generated  three  cases  of  Abreast, 
Single-file  and  Staggered  target  formations.  In  all  these  cases  the  radius  of  the  array  was  4  ft.  For  the 
Abreast  case,  the  two  targets  are  originally  at  (-300m, -300m)  and  (-300m, -350m).  For  the  Single-file  case, 
the  two  targets  are  originally  at  (-300m, 100m)  and  (-250m, 100m).  In  Staggered  formation  there  are  three 
targets  at  initial  positions  (-300m, 100m),  (-250m, 100m),  (-275m, 150m).  In  all  these  cases,  the  targets  move 
to  East  with  a  constant  speed  of  25m/s.  The  signal-to-noise  ratio  (SNR)  was  lOdb  for  these  cases.  The 
results  for  these  three  cases  are  described  in  details  below. 

•  Single-file  case:  The  ground  truth  plot  associated  with  the  two  targets  is  shown  in  Figure  33. 

The  set  of  the  principal  components  (PC’s)  based  on  the  eigenvalue  decomposition  algorithm  for  the 
’’window”  index  5  are  also  shown  in  this  figure.  In  all  the  figures  shown  below  the  left  column  depicts 
the  raw  data  from  the  five  sensors,  the  middle  column  depicts  the  centered  data  as  discussed  in  the 
preprocessing  Section  6.4.1,  and  the  right  column  represents  the  independent  components  that  are 
estimated  by  the  Fast  ICA  algorithm.  The  ICA  is  performed  in  two  stages  where  in  the  first  stage 
the  algorithm  is  allowed  to  estimate  all  the  independent  components  available  in  the  data  set.  In 
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the  second  stage  the  algorithm  is  imposed  to  estimate  two  or  three  independent  components  based 
on  the  dominant  PC  A  singular  values  as  shown  in  Figure  33.  Figure  34  shows  the  ICA  of  the  data 
set  for  the  window  index  5  with  estimation  of  five  components.  It  can  be  observed  from  this  figure 
that  there  are  clearly  two  source  signals  and  the  other  three  are  noise  sources.  In  the  next  set  of 
simulations  instead  of  estimating  five  components,  the  Fast  ICA  algorithm  is  asked  to  generate  only 
two  components.  These  results  are  also  shown  in  Figure  34  due  to  space  limitations.  Specifically, 
the  1st  and  the  2nd  signals  from  the  top  in  the  right  column  is  the  response  of  the  Fast  ICA  that  is 
limited  to  generating  two  independent  components.  This  should  confirm  clearly  the  presence  of  only 
two  targets. 

•  Abreast  case:  The  ground  truth  plot  associated  with  the  two  targets  is  shown  in  Figure  35.  The 
set  of  the  principal  components  based  on  the  eigenvalue  decomposition  algorithm  for  the  the  window 
index  14  are  also  shown  in  this  figure.  Figure  36  shows  the  ICA  of  the  data  set  for  the  window  index 
14  with  estimation  of  five  components.  It  can  be  observed  from  these  figures  that  there  are  clearly 
two  source  signals  and  the  other  three  are  noise  sources.  In  the  next  set  of  simulations  instead  of 
estimating  five  components  the  Fast  ICA  algorithm  is  asked  to  generate  only  two  components.  These 
results  are  also  shown  in  Figure  36  (the  4th  and  the  5th  independent  signals  in  the  right  column) 
which  should  again  confirm  the  presence  of  only  two  targets. 


Figure  33:  Ground  truth  and  Principal  components  (PCA)  based  on  Eigenvalue  Decomposition  (EVD) 
algorithm  for  the  Single-File  formation  case  and  for  the  window  index  5. 


•  Staggered  case:  The  above  results  are  now  repeated  for  the  Staggered  synthetic  data.  The 

ground  truth  plot  associated  with  the  three  targets  is  shown  in  Figure  37.  The  set  of  the  principal 
components  based  on  the  eigenvalue  decomposition  algorithm  for  the  window  index  40  are  also  shown 
in  this  figure.  Figures  38  show  the  ICA  of  the  data  set  for  the  window  index  40  with  estimation  of 
five  components.  It  can  be  observed  from  this  figure  that  there  are  clearly  three  source  signals  and 
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Figure  34:  Fast  ICA  algorithm  applied  to  Single-File  formation  data.  The  right  column  indicates  the  inde¬ 
pendent  components  using  five  components  for  window  index  5.  Limiting  the  algorithm  to  two  independent 
components  results  in  the  1st  and  2nd  signals  from  the  top  right  column  as  the  independent  components. 
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Figure  35:  Ground  truth  and  Principal  components  (PCA)  based  on  Eigenvalue  Decomposition  (EVD) 
algorithm  for  the  Abreast  formation  case  and  the  window  index  14. 


the  other  two  are  noise  sources.  In  the  next  set  of  simulations  instead  of  estimating  five  components 
the  Fast  ICA  algorithm  is  asked  to  generate  only  three  components.  These  results  are  also  shown  in 
Figure  38  (the  2nd ,  3rd  and  the  5th  independent  components  in  the  right  column)  which  should  again 
confirm  the  presence  of  only  three  targets. 

6.6  Simulation  Results  of  Fast  ICA  on  SAFE  II  Data 

The  ICA  is  performed  for  each  window  of  1  second  duration  and  for  each  three  (3)  nodes  consisting  of  five 
(5)  sensors.  The  first  data  set  (503  series)  chosen  here  has  two  independent  sources.  The  ground  truth  plots 
associated  with  the  three  nodes  are  shown  in  Figure  39.  As  in  the  previous  set  of  results,  the  horizontal 
axes  in  this  figure  represents  the  55  window”  index.  The  fast  ICA  is  performed  on  this  data  set  in  two  stages 
where  in  the  first  stage  the  algorithm  is  allowed  to  estimate  all  the  independent  components  available 
in  the  data  set.  In  the  second  stage  the  algorithm  is  imposed  to  estimate  the  independent  components 
based  on  the  dominant  PCA  singular  values  (see  Figure  40).  In  all  the  figures  shown  below  related  to  the 
Fast  ICA  independent  components  the  left  column  depicts  the  raw  data  from  the  five  sensors,  the  middle 
column  depicts  the  centered  data  as  discussed  in  the  preprocessing  Section  6.4.1,  and  the  right  column 
represents  the  independent  components  that  are  estimated.  Figure  41  shows  the  ICA  of  data  set  for  node 
1  with  estimation  of  five  components,  where  our  goal  was  to  initially  estimate  five  components.  From  this 
figure  the  presence  of  two  source  signals  and  three  noisy  signals  is  rather  evident.  To  substantiate  this,  the 
Fast  ICA  algorithm  is  run  by  requiring  it  to  generate  only  two  independent  components  (source  signals). 
The  results  are  shown  directly  in  Figure  41  (due  to  space  limitations)  corresponding  to  the  1st  and  the  3rd 
signals  from  the  top  in  the  right  column  subfigure,  which  demonstrates  the  independence  of  the  two  sources 
for  this  case.  Very  similar  results  for  nodes  2  and  3  are  also  obtained,  however  due  to  space  limitations  are 
not  shown  here. 
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Figure  36:  Fast  ICA  algorithm  applied  to  Abreast  formation  data.  The  right  column  indicates  the  indepen¬ 
dent  components  using  five  components  for  window  index  14.  Limiting  the  algorithm  to  two  independent 
components  results  in  the  4th  and  5th  signals  from  the  top  right  column  as  the  independent  components. 
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Figure  37:  Ground  truth  and  Principal  components  (PCA)  based  on  Eigenvalue  Decomposition  (EVD) 
algorithm  for  the  Staggered  formation  case  and  the  window  index  40. 


Next,  in  Figures  42  we  have  repeated  the  above  experiments  by  using  the  window  index  130  for  nodes 
2  &  3  and  index  90  for  node  1  (refer  to  Figure  39  for  detailed  references).  As  can  be  seen  from  Figure  39 
the  above  window  indexes  correspond  to  the  situation  where  the  moving  target  is  at  its  closest  proximity 
to  the  sensors.  We  are  particularly  interested  in  this  situation  in  order  to  determine  the  performance  of 
the  Fast  ICA  algorithm  with  respect  to  its  dependence  to  the  distance  of  targets  to  sensors.  It  can  be 
shown  from  Figure  42  that  the  dominance  of  two  targets  is  more  profoundly  indicated  for  data  newr2503 
and  newr3503  as  compared  to  the  results  that  we  obtained  in  Figure  40  for  data  newr2503  to  newr3503 
for  window  indexes  20  and  30,  respectively.  These  results  translate  further  into  better  resolution  of  source 
signals  from  noisy  signals  as  shown  in  Figure  43.  The  results  using  two  independent  sources  demonstrate 
(the  1st  and  the  3rd  signals  from  the  top  right  column)  the  fact  that  there  are  indeed  two  source  signals  in 
this  case.  Very  similar  results  are  also  obtained  for  node  2  but  due  to  space  limitations  are  not  included 
here.  However,  the  results  for  data  newrl503  (node  1)  for  the  window  index  130  is  inaccurate,  as  the  Fast 
ICA  algorithm  is  suggesting  the  presence  of  five  independent  sources.  This  unexpected  behavior  is  an  issue 
that  we  need  to  further  investigate  in  Phase  II. 

The  next  data  set  chosen  is  from  the  510  series  where  one  distinctive  source  signal  is  present.  Figure 
44  depicts  the  ground  truth  plots  corresponding  to  nodes  1  to  3.  Based  on  the  PCAs  obtained  associated 
with  the  window  index  20  for  all  three  nodes,  one  can  infer  for  node  1  the  presence  of  two  sources,  whereas 
nodes  2  3  confirm  the  presence  of  one  source  signal.  To  see  this  Figure  45  shows  the  separated  signals 

for  node  2  using  five  components.  This  figure  shows  components  that  are  very  similar  to  those  obtained 
in  node  3  (not  shown  due  to  space  limitations).  A  careful  analysis  of  the  results  for  node  1  does  suggest 
consistency  with  PCA  results.  In  other  words,  there  are  possibly  two  independent  source  signals.  However, 
this  conclusion  is  not  supported  by  the  results  seen  from  Figure  45  and  those  obtained  from  node  3  where 
there  is  an  indication  of  only  one  independent  source  signal.  Limiting  the  Fast  ICA  to  generate  only  one 
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Figure  38:  Fast  ICA  algorithm  applied  to  Staggered  formation  data.  The  right  column  indicates  the 
independent  components  using  five  components  for  window  index  40.  Limiting  the  algorithm  to  three 
independent  components  results  in  the  2nd,  3rd  and  the  5th  signals  from  the  top  right  column  as  the 
independent  components. 
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Figure  40:  Principal  components  (PCA)  based  on  Eigenvalue  Decomposition  (EVD)  algorithm  applied  to 
SAFE  II  real  data  newr!503  to  newr3503  for  the  window  indexes  20,  20,  and  30,  respectively. 


independent  component  results  in  the  4th  signal  from  the  top  in  the  right  column  of  Figure  45  as  being  the 
source  signal. 

The  final  data  set  chosen  is  from  the  537  series  where  five  distinctive  source  signals  are  present.  Results 
in  Figure  46  depicts  the  ground  truth  associated  with  nodes  1  to  3  and  the  definitions  corresponding  to 
the  window  indexes.  Figure  47  depicts  the  PCA  estimates  for  window  index  125  for  nodes  2  &  3  and  index 
100  for  node  1.  Although  node  3  clearly  indicates  the  presence  of  five  dominant  signals,  this  is  to  a  lesser 
extent  confirmed  by  node  2  and  even  to  a  lesser  degree  by  node  1.  Figures  48  depicts  the  five  independent 
components  that  are  generated  by  the  Fast  ICA  algorithm  for  node  2,  which  does  suggest  the  presence  of 
only  two  source  signals  as  opposed  to  five  components.  The  reasons  for  this  discrepancy  and  error  will  be 
investigated  in  more  details  in  Phase  II. 

The  next  data  set  chosen  from  the  537  series  are  from  the  window  index  200  for  node  1  and  index  250 
for  nodes  2  &  3.  By  plotting  the  PCA  values  in  this  case  it  can  be  seen  (deleted  due  to  space  limitations) 
that  nodes  1  &  3  clearly  indicate  the  presence  of  five  dominant  signals,  however  this  is  to  a  lesser  extent 
confirmed  by  node  2.  Figures  49  depicts  the  five  independent  components  that  are  generated  by  the  Fast 
ICA  algorithm  for  node  2.  Finally,  from  the  result  of  the  PCA  estimates  for  the  window  index  400  for 
all  the  three  nodes,  it  follows  that  nodes  2  &  3  indicate  the  presence  of  five  independent  source  signals, 
whereas  node  1  does  only  suggest  the  presence  of  one  independent  source  signal.  The  corresponding  Fast 
ICA  generated  independent  components  for  node  1  are  shown  in  Figure  50,  again  indicating  an  error  in 
source  separation.  In  contrast  the  results  in  nodes  2  &  3  both  indicate  the  presence  of  5  source  signals, 
but  these  plots  are  omitted  due  to  space  limitations.  1 
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Figure  41:  Fast  ICA  algorithm  applied  to  SAFE  II  real  data  newrl503.  The  right  column  indicates 
the  independent  components  using  five  components  for  window  index  20.  Limiting  the  algorithm  to  two 
independent  components  results  in  the  1st  and  3rd  signals  from  the  top  right  column  as  the  independent 
components. 
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Figure  42:  Principal  components  (PCA)  based  on  Eigenvalue  Decomposition  (EVD)  algorithm  applied  to 
SAFE  II  real  data  newr!503  to  newr3503  for  the  window  indexes  130,  130,  and  90,  respectively. 
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Figure  43:  Fast  ICA  algorithm  applied  to  SAFE  II  real  data  newr2503.  The  right  column  indicates  the 
independent  components  using  five  components  for  window  index  130.  Limiting  the  algorithm  to  two 
independent  components  results  in  the  1st  and  3rd  signals  from  the  top  right  column  as  the  independent 
components. 
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Figure  44:  Ground  truth  diagrams  for  SAFE  II  real  data  newr!510  to  newr3510. 
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Figure  45:  Fast  ICA  algorithm  applied  to  SAFE  II  real  data  newr2510.  The  right  column  indicates 
the  independent  components  using  five  components  for  window  index  20.  Limiting  the  algorithm  to  one 
independent  component  results  in  the  4th  signal  from  the  top  right  column  as  the  independent  component. 
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Figure  46:  Ground  truth  diagrams  for  SAFE  II  real  data  newr!537  to  newr3537 


Figure  47:  Principal  components  (PCA)  based  on  Eigenvalue  Decomposition  (EVD)  algorithm  applied 
SAFE  II  real  data  newr!537  to  newr3537  for  window  indexes  100,  125,  and  125,  respectively. 
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Figure  48:  Fast  ICA  algorithm  applied  to  SAFE  II  real  data  newr2537.  The  right  column  indicates  the 
independent  components  using  five  components  for  window  index  125. 
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Observations: 


•  The  PCA  analysis  conducted  shows  that  the  degree  of  dominant  features  could  be  different  depending 
on  the  particular  node  and  the  location  of  the  targets  to  sensors  that  one  considers.  Specifically,  for 
the  510  series  data  it  can  be  argued  that  the  sets  newr2510  and  newr3510  may  be  approximated 
by  one  independent  component  whereas  the  newrl510  data  suggested  two  strongly  dominant  source 
signals.  This  is  interesting,  as  we  clearly  know  that  there  is  only  one  target  in  this  case,  therefore 
it  would  be  interesting  to  determine  the  nature  of  the  second  feature  present  in  this  node,  as  well 
as  to  use  a  data  fusion  algorithm  to  integrate  these  results  to  generate  valid,  verifiable,  and  reliable 
conclusions. 

•  The  above  observation  also  applies  to  series  537  where  in  nodes  1  &  2  in  window  indexes  100  and 
125  there  seem  to  exist  less  than  five  independent  source  signals.  This  is  interesting  as  we  know 
that  there  are  five  independent  source  signals.  Integration  and  fusion  of  these  specific  nodes  with 
the  other  nodes  should  be  able  to  remove  this  ambiguity  and  inconsistency.  This  is  a  topic  of  future 
research  and  investigation  in  Phase  II. 

•  In  all  the  simulation  results  conducted  above  it  is  interesting  to  note  that  the  uniformly  “flat” 
component  signal  is  the  most  dominant  pattern  in  all  the  cases  and  for  all  the  nodes.  We  are 
speculating  (without  having  any  further  information)  that  this  should  correspond  to  a  “feature”  that 
is  very  common  for  all  the  vehicles  present  in  real  data. 

7  Composite  Beampattern 

In  practice,  the  arrays  are  normally  deployed  hundreds  of  meters  apart.  However,  the  farther  the  arrays 
are,  the  less  coherent  the  signals  received  from  the  arrays.  As  we  know,  most  of  DOA  detectors  rely  on 
coherence  of  the  signals  received  from  different  sensors.  As  a  consequence,  the  original  version  of  the 
ESPRIT  algorithm  described  in  the  previous  section  that  uses  the  signals  from  all  the  arrays  in  the  DOA 
estimation  process  will  not  directly  apply  to  this  case. 

To  deal  with  this  problem,  an  alternative  method  that  relies  on  forming  the  “composite  beampattern” 
is  proposed  here.  The  motivation  behind  this  method  is  that  even  though  one  array  may  not  resolve  two 
closely  spaced  targets,  the  two  targets  might  look  different  in  different  arrays.  Additionally,  even  if  one 
array  cannot  detect  all  the  targets,  the  other  arrays  will  pick  up  these  missed  detected  targets.  For  example, 
assume  that  the  arrays  are  linear  and  the  two  targets  form  a  vertical  line  with  respect  to  one  array.  Now, 
if  the  second  array  is  placed  perpendicular  to  the  previous  array,  it  will  see  a  parallel  line  formed  by  the 
targets.  This  property,  in  turn,  will  result  in  a  slight  difference  in  the  peaks  of  the  beampatterns  formed 
by  different  arrays.  Thus,  the  composite  beampattern  of  all  the  arrays  would  exhibit  peaks  corresponding 
to  all  the  detected  targets.  The  proposed  procedure  is  as  follows: 

Let  &i(0),  &2(0)  bs(0)  and  b^[6)  be  the  individual  beampatterns  of  arrays  1  2,3  and  4,  respectively. 
Suppose  that  for  a  target  9\ ,02,03  and  #4  are  the  DOA’s  relative  to  the  center  of  each  array,  as  shown  in 
Figure  5.  Let  us  denote  the  composite  beampattern  by  bx(a)  where  a  is  the  DOA  with  respective  to  the 
origin  and  take  the  DOA,  0i,  detected  by  array  1  as  the  reference.  The  known  position  of  the  center  of 
array  1  is  assumed  to  be  at  {x\7yi).  Now,  if  the  target  has  a  DOA  of  9\  relative  to  the  center  of  array  1 
and  a  DOA  of  a  relative  to  the  origin,  by  knowing  the  position  of  array  1  (denoted  by  (xifyi)),  we  could 
represent  the  position  of  the  target,  (u,v)7  as, 

yi  —  x\  tan  9\ 

u  —  - — 

tan  a  —  tan  9\ 

v  =  u  tan  a 
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Figure  50:  Fast  ICA  algorithm  applied  to  SAFE  II  real  data  newr!537  using  five  components. 
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for  84  <  a  <  tan-1  |j-.  After  (u,v)  is  obtained,  one  can  easily  represent  the  corresponding  DOA’s  82,  83 
and  84,  relative  to  the  center  of  arrays  2,3  and  4,  respectively,  i.e. 

8i  —  tan-1  — — —  *  =  2,3,4 

U  —  Xi 

Then,  the  composite  beampattern  at  DOA  a  is  bx(a)  =  82(82) +  ki(8^)  + 84(84).  Note  that  since  we  used 
beampattern  1  as  the  reference,  84  is  actually  fixed  (i.e.  the  peak  of  beampattern  1)  and  hence  84(84)  is  also 
fixed.  Thus,  it  does  not  make  any  difference  whether  or  not  we  add  bl(8i)  to  the  composite  beampattern. 

The  procedure  involves  incrementing  a  within  the  range  84  <  a  <  tan-1  and  forming  the  composite 
beampattern.  The  peaks  of  the  composite  beampattern  correspond  to  the  DOA’s  of  all  the  targets.  It 
must  be  pointed  out  that  a  is  incremented  with  the  same  resolution  as  that  used  to  form  the  individual 
beampattern  for  each  array.  This  due  to  the  fact  that  the  resolution  of  the  composite  beampattern  is 
determined  by  the  individual  resolutions. 


Figure  51:  Forming  Composite  Beampattern. 


One  of  the  main  goals  to  pursue  for  the  next  phase  of  this  research  effort  is  to  experimentally  implement 
and  test  the  above  algorithms  for  both  a  real  (e.g.  SAFE  II)  data  set  as  well  as  our  synthetically  generated 
data  sets  consisting  of  several  closely  spaced  targets  in  single-file,  staggered  and  abreast  formations.  In 
particular,  our  Phase  II  efforts  should  address  the  following  problems  of  the  above-mentioned  algorithm. 

•  The  spatial  invariance  approach  using  ESPRIT  relies  on  coherence  of  the  signals  arriving  at  the 
sparse  arrays.  However,  this  assumption  is  not  valid  when  the  spacing  A  is  very  large  or  in  presence 
of  wind  noise  and  cluttered  environment.  Thus,  this  method  needs  to  be  extended  to  address  the 
very  issue  in  our  specific  application. 
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•  Presently,  the  spatial  invariance  approach  using  ESPRIT  is  applied  to  narrowband  sources.  We  plan 
to  develop  a  subband  version  of  this  algorithm  for  wideband  sources, 

•  Implementation,  testing  and  bench  marking  of  the  composite  beampattern  algorithm  on  both  real 
and  synthetic  data,  and 

•  Comparison  of  the  performance,  robustness  and  accuracy  of  the  approaches  presented  above  for 
several  spatially  invariant  subarray  configurations. 


8  Conclusions  &;  Future  Plans  for  Phase  II 

The  problem  of  detection,  localization  and  tracking  of  multiple  closely  spaced  targets  from  several  arrays 
of  unattended  acoustic  sensors  was  considered  in  this  Phase  I  project.  The  research  activities  in  this  phase 
were  concentrated  on  the  development  of  various  high-resolution  and  robust  algorithms  for  detecting  and 
resolving  multiple  targets  in  tight  formations.  Barrier  issues  critical  to  the  development  of  a  realistic 
acoustic  signature  analysis  system  were  also  investigated.  Several  wideband  DOA  estimation  methods, 
namely  incoherent  wideband  MUSIC  (IWM),  coherent  steered  covariance  matrix  (STCM),  and  weighted 
subspace  fitting  (WSF)  methods,  were  carefully  studied  (Section  4)  and  benchmarked  on  both  synthesized 
as  well  as  some  limited  real  data  (SAFE  II).  Our  studies  indicated  that  the  STCM  is  an  efficient  method 
for  locating  and  tracking  target  formations  though  it  may  not  be  very  effective  in  resolving  the  individual 
targets,  especially  when  they  are  very  close  together.  On  the  other  hand,  WSF  method  exhibited  very 
good  resolution  in  separating  different  targets  in  a  formation.  The  results  on  abreast,  staggered  and 
single-file  formations  revealed  the  excellent  resolution  of  this  method.  However,  this  algorithm  requires  an 
exhaustive  and  laborious  search  though  the  entire  space  in  order  to  reduce  the  ambiguity  in  estimating  the 
DOA’s.  Our  simulation  results  also  demonstrated  an  interesting  observation  that  if  STCM  is  first  applied 
as  a  preprocessor  to  detect  and  locate  the  target  formations,  then  WSF  could  efficiently  be  implemented 
on  a  reduced  search  space  to  resolve  the  individual  targets.  This  is  one  of  the  items  that  we  intend  to 
research  further  in  Phase  II.  Another  methodology  for  separating  the  targets  (sources)  was  also  extensively 
studied.  This  approach  is  based  upon  blind  source  separation  using  ICA  methods.  Two  different  ICA 
algorithms  namely  IICA  and  Fast  ICA  were  investigated  and  applied  to  this  problem.  Our  results  revealed 
that  from  both  computational  simplicity  and  performance  Fast  ICA  is  a  preferred  algorithm  for  source 
separation.  This  algorithm  is  very  simple  and  easy  to  use  comparing  to  other  ICA  methods.  One  of  the 
objectives  in  Phase  II  will  be  to  investigate  and  develop  ”  nonlinear  version  of  the  Fast  ICA  and  investigate 
its  performance  as  compared  to  the  linear  ICA  algorithms  investigated  in  Phase  I.  The  sensitivity  of  both 
the  linear  ICA  and  the  nonlinear  ICA  algorithms  to  the  variations  in  the  nature  of  the  sources,  choice  of 
the  contrast  functions  employed,  noise  and  other  complicating  attributes  will  be  fully  investigated  in  Phase 
II. 


Additionally,  we  investigated  the  applicability  of  the  spatial  invariance  scheme  using  a  rectangular  array 
consisting  of  subarrays /elements  of  the  nodes  that  are  sparsely  distributed  in  the  field.  The  expected  benefit 
of  this  approach  is  wider  aperture  and  hence  better  resolution.  Preliminary  results  of  this  narrowband 
algorithm  demonstrated  the  potential  of  the  method  for  our  problem.  Finally,  we  proposed  a  new  algorithm 
for  combining  the  estimated  DOA’s  from  different  arrays  in  order  to  yield  unambiguous  target  detection 
and  tracking. 

8.1  Future  Plans  for  Phase  II 

We  shall  continue  our  current  research  activities  in  Phase  I  toward  developing  dedicated  high-resolution 
and  robust  DOA  estimation  methods  for  detection,  separation  and  subsequent  classification  and  tracking 
of  multiple  closely  spaced  targets.  Various  barrier  issues  and  major  impediments  to  the  development  of 
practical  and  real-life  systems  have  been  identified  in  Phase  I  of  this  project.  Among  the  major  issues  that 
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we  have  identified  for  further  investigation  and  research  are:  (a)  solutions  to  the  above-mentioned  problems 
under  non-uniform  and  arbitrary  geometries,  taking  into  account  the  advantages  and  disadvantages  of 
multiple  invariance  geometries  such  as  the  planar  geometry  that  is  currently  being  utilized  by  the  US 
Army  Labs,  (b)  robustness  and  sensitivity  analysis  of  the  proposed  as  well  as  standard  methods  to  the 
correlated  and  non-Gaussian  noise  sources,  wind,  heat/cold  and  environmental  effects  and  competing 
clutter,  as  well  as  possible  losses  in  data  and  failure  in  the  array  elements  and/or  array  nodes  due  to 
unexpected  circumstances,  (c)  multiple  target  tracking  and  classification  for  unattended  ground  sensors, 
and  (d)  ICA  and  blind  source  separation  techniques  for  determining  the  number  of  sources  from  their 
convolutive  mixture.  In  particular,  Phase  II  research  involves  development  of  algorithms  that  address 
all  the  above-mentioned  requirements  and  can  be  integrated  as  part  of  the  Government  acoustic  data 
acquisition  and  analysis  system.  The  testing,  performance  evaluation  and  benchmarking  of  the  algorithms 
will  be  performed  on  a  large  number  of  multiple  closely  spaced  target  scenarios  under  various  realistic 
operating  conditions.  The  specific  research  objectives  of  the  Phase  II  research  study  are  listed  below. 

1.  Considering  the  geometry  of  the  arrays  and  the  nodes  (typically  nonuniform)  commonly  used  in  bat¬ 
tlefield  scenarios,  our  goal  is  to  develop  efficient  dedicated  DOA  estimation,  source  separation  and 
classification  algorithms  for  high-resolution  and  robust  target  bearing  angle  estimation  and  separa¬ 
tion.  The  purpose  here  is  to  develop  novel  ideas  for  optimally  combining  different  elements  of  the 
nodes  in  order  to  substantially  improve  spatial  resolution  for  distinguishing  targets  in  abreast,  stag¬ 
gered  or  single-file  formations.  The  proposed  methods  use  spatial  invariance  by  taking  a  rectangular 
array  consisting  of  sub  arrays /elements  of  the  nodes  that  are  sparsely  distributed  in  the  field.  The 
benefit  is  certainly  wider  aperture  and  much  finer  resolution.  To  avoid  cyclical  ambiguity  due  to 
spatial  aliasing,  we  propose  a  new  way  of  combining  the  sparse  arrays  and  using  subband  processing. 
Furthermore,  by  developing  MUSIC  and  ESPRIT-based  approaches  for  this  architecture  it  is  our  goal 
to  improve  on  these  standard  techniques.  We  are  also  interested  in  considering  the  advantages  and 
disadvantages  of  various  wideband  DOA  estimation  techniques,  such  as  STCM-based  method,  and 
WSF  technique,  to  name  a  few.  The  issues  to  address  here  will  deal  with  the  design,  development  and 
simulation  of  the  proposed  techniques  suitable  for  the  specific  geometry(ies)  that  are  used  in  real-life 
battlefield  scenarios.  The  primary  criterion  is  to  develop  algorithms  that  can  be  implemented  in 
real-time.  In  addition,  the  methods  should  not  rely  on  prior  knowledge  of  the  number  of  targets, 
target’s  dynamical  information  and  initial  conditions,  and  background  interference  and  clutter. 

2.  We  will  continue  our  present  work  on  subband  MUSIC  algorithm  and  extend  the  methodologies  for 
different  array  combinations  for  multiple  node  configurations.  This  method  can  also  be  combined  with 
the  spatial  invariance  methods  as  described  above.  The  subband  processing  and  DOA  estimation  can 
substantially  reduce  the  ambiguity  due  to  the  spatial  aliasing  of  sparse  arrays  while  giving  accurate 
DOA’s.  We  shall  investigate  how  to  resolve,  with  minimum  ambiguity  and  high  confidence,  the 
DOA’s  associated  with  multiple  closely  spaced  targets  in  spatially  aliased  data  due  to  sparse  nature 
of  the  arrays  employed  for  this  problem.  This  is  especially  important  in  presence  of  spatially  colored 
interference  and  other  interferences  such  as  wind  noise  and  environmental  factors. 

3.  We  will  continue  on  our  preliminary  study  on  using  ICA  as  a  tool  to  perform  blind  source  separation. 
The  goal  here  is  to  investigate  whether  or  not  in  tight  target  formations  ICA  could  potentially 
determine  the  number  of  sources.  Although,  this  method  does  not  have  the  tracking  ability,  it  is 
simple  and  does  not  suffer  from  the  resolution  versus  aliasing  tradeoff  as  with  the  DOA  estimation 
methods.  The  practical  suitability  for  this  technique  to  our  particular  problem  becomes  indeed  quite 
challenging  when  correlated  noise  and/or  multipath  is  present. 

4.  We  would  also  investigate  the  issue  of  track  association  and  target  classification.  Track/data  associ¬ 
ation  can  be  used  to  remove  the  false  peaks  in  the  DOA  estimates  and  obtain  smooth  tracks.  This  is 
particularly  important  in  staggered  formations  when  the  tracks  may  be  difficult  to  associate.  Once 
this  is  accomplished,  based  upon  the  tonal/movement  features  of  the  targets  classification  can  be 
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made.  The  recursive  high  order  correlation  method  and  sequential  Bayes  classifier  developed  by  the 
PI  will  be  employed  for  track  association  and  target  classification  tasks. 

5.  The  algorithms  will  be  tested  and  analyzed  on  new  data  sets  that  will  be  acquired  from  the  US  Army 
TACOM-ARDEC.  These  data  sets,  which  will  be  collected  in  different  environmental  conditions, 
include  several  multiple  target  cases  that  involve  various  target  scenarios  and  high  density  of  clutter 
and  correlated  interference.  Additionally,  we  continue  to  generate  several  sets  of  realistic  synthesized 
data  by  inserting  multiple  closely  spaced  targets  in  different  tight  formations.  This  would  enable  us 
to  study  the  performance  of  the  overall  system  in  terms  of  accuracy,  robustness  and  stability  of  the 
DOA  estimation  algorithms.  Moreover,  we  shall  continue  to  investigate  how  the  algorithms  perform 
under  missing  data  or  sensor  failure  conditions.  This  is  extremely  critical  in  real  battlefield  situations 
and  in  adverse  environmental  conditions.  We  would  like  to  see  how  our  physical  understanding  of 
the  behavior  and  properties  of  the  array-node  configurations  can  inform  the  choice  of  algorithms  to 
remedy  these  difficult  problems. 
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